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1.  Summary 

The  main  purpose  of  the  Phase  I  effort  was  to  set  up  the  general  framework  for  the  wavefront  (  WF) 
evolution  technique  and  its  initial  implementation  The  framework  we  developed  includes  algorithms 
for  evolution  of  WFs  and  fields  associated  with  them,  as  well  as  algorithms  for  computation  of  the 
scattered  fields  (either  directly,  or  through  the  prior  computation  of  surface  currents).  An  emphasis 
was  put  on  incorporating,  in  a  fashion  consistent  with  reciprocity,  edge  diffraction  effects  into  the 
WF  approach.  A  UTD  formulation,  modified  to  reproduce  rigorously  the  canonical  solution  in  the 
near  and  transition  regions,  has  been  constructed  [1]  to  describe  the  “non-optical  character  of 

diffracted  rays,  and  to  evaluate  diffracted  fields. 

An  additional  benefit  of  our  development  of  the  WF  method  is  that  it  will  provide  an  essential 
element,  a  general  parametric  representation  of  the  Ansatz ,  to  the  asymptotic  high  frequency  integral 
equation  formulations  [2, 3]  based  on  the  stationary  point  method.  In  other  words,  the  HF  code  we 
plan  to  construct  may,  as  a  byproduct,  provide  an  input,  for  a  general  geometry ,  to  an  asymptotic 
high-frequency  integral  equation  solver. 

Phase  I  results  may  form  the  basis  for  the  further  development  of  an  efficient  high-frequency 
scattering  code. 


2.  Development  of  the  WF  evolution  framework 

2.1.  Main  features  of  WF  propagation  approach 

In  the  context  of  high-frequency  asymptotic  expansions  the  WF  (or  a  phasefront)  is  defined  math¬ 
ematically  as  a  surface  on  which  the  considered  propagating  fields  have  a  constant  phase.  The  set 
of  all  WFs  forms  a  one-parameter  family  of  surfaces,  characterized  by  the  phase  or,  equivalently, 
by  the  length  of  the  path  measured  from  the  wave  source.  In  the  following  we  use  the  path  length 
as  the  evolution  parameter  of  the  WF . 

Most  of  the  recent  activity  in  the  development  WF  evolution  techniques  occurred  in  the  area 
of  geophysical  (seismic)  problems  [4, 5, 6, 7, 8].  In  such  application  the  main  interest  lies  in  describ¬ 
ing  wave  propagation  through  penetrable  inhomogeneous  media,  i.e.,  in  the  wave  refraction.  In 
electromagnetics,  on  the  other  hand,  the  more  relevant  phenomena  are  reflection  and  diffraction. 
Correspondingly,  our  WF  evolution  algorithms  concentrate  primarily  on  these  processes. 

Wavefront  evolution  may  be  considered  as  an  improved  ray-tracing  method,  in  which  the  number 
of  rays  is  dynamically  adjusted ,  while  the  WF  expands  or  shrinks.  New  rays  are  also  created  at 
shadow  and  reflection  boundaries,  and  in  diffraction  processes. 

The  advantage  of  the  WF  evolution  approach,  compared  to  the  conventional  ray  tracing,  is 
that  it  maintains  an  approximately  constant  ray-ray  spacing  h ,  even  in  processes  involving  multiple 
reflection  and  diffraction;  this  is  achieved  due  to  the  dynamic  adjustment  of  the  number  of  rays. 
In  our  implementation  the  WF  is  represented  as  a  meshed  (triangulated)  surface  with  well  defined 
boundaries;  consequently,  adjustment  of  the  number  of  rays  is  realized  as  mesh  simplification  or 
complexification.  The  new  rays  maintain  a  precise  WF  definition,  due  to  the  use  of  the  local 
second-order  curvature-based  surface  representations. 
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Another  advantage  of  our  implementation  of  the  WF  evolution  algorithms  is  that  they  allow 
a  precise,  second  order  in  h  determination  of  shadow  and  reflection  boundaries',  this  is  also  achieved 
by  creating  additional  rays  located  precisely  at  these  boundaries. 

Probably  the  most  significant  benefit  of  the  WF  evolution  methods  is  that  they  allow,  in  a 
natural  way,  for  the  consistent  inclusion  of  multiple  diffraction  mechanisms  and  diffraction  occurring 
in  combination  with  reflection  processes.  As  we  describe  below  in  more  detail,  in  our  Phase  I  work 
we  initiated  the  implementation  of  the  geometrical  construction  of  diffracted  WFs  due  to  edge 
diffraction,  and  the  corresponding  evaluation,  according  to  the  UTD  asymptotic  expansion,  of  the 
fields  associated  with  the  diffracted  WFs. 

Our  WF  evolution  algorithms  are  currently  implemented  for  scatterers  defined  by  some  se¬ 
lected  analytically  defined  surfaces,  or  by  facetized  (triangulated)  surfaces.  They  can  be,  however, 
extended  in  a  relatively  straightforward  way  to  general  parametric  surfaces,  defined,  e.g.,  in  terms 
of  bicubic  patches.  The  only  new  element  required  for  such  a  generalization  is  a  set  of  func¬ 
tions  providing  information  on  the  normals  and  curvatures  of  the  scatterer  surface,  and  computing 
intersections  of  lines  with  the  surface. 

We  note  that  our  developments  are  based  on  the  “Lagrangian”  formulation  of  the  WF  the¬ 
ory,  in  which  WF  are  characterized  by  geometrical  and  field  parameters  associated  with  a  set  of 
vertices  (“markers”)  moving  in  space  with  the  changing  evolution  parameter.  There  have  been 
recent  developments  in  the  alternative  “Eufernn”  methods  (Refs.  [9, 10, 11]  and  private  communi¬ 
cations),  in  which  the  WF  parameters  are  associated  with  points  of  a  fixed  spatial  grid  (possibly  in 
the  higher  number  of  dimensions).  These  methods  have  certain  advantages  in  terms  of  simplicity 
of  formulation,  and  provide  a  natural  way  of  defining  WF  with  an  approximately  constant  spatial 
resolution  (equal  to  the  Eulerian  grid  spacing).  The  latter  feature,  however,  may  become  a  hin¬ 
drance  in  describing  the  WF  behavior  near  caustics,  because  the  Eulerian  formulation  naturally 
provides  only  storage  for  a  limited  amount  of  information  per  unit  volume.  Therefore,  situations 
where  rays  converge  require  special  treatments,  such  as  using  higher  dimensions  or  extension  in 
the  evolution  parameter.  We  note  that  the  WF  behavior  near  caustics  is  also  relevant  in  describing 
diffraction,  since  diffraction  sources  are  always  caustics  of  the  diffracted  WF. 

It  should  be  noted  that  neither  Lagrangian  nor  Eulerian  formulations  provide,  per  se,  any  new 
physical  mechanism  or  any  new  mathematical  apparatus  for  describing  GO  evolution  or  diffrac¬ 
tion  of  waves.  These  processes  have  to  modeled,  in  both  approaches,  using  the  asymptotic  high 
frequency  methods.  Also,  geometrical  problems  associated  with  the  field  propagation,  such  as 
detection  of  a  scatterer  or  identification  diffraction  sources  on  the  scatterer,  have  to  be  resolved 
independently  of  the  use  of  a  Lagrangian  or  an  Eulerian  approach. 

For  all  the  reasons  listed  above  we  decided  to  stay,  in  our  developments,  with  the  more  con¬ 
ventional  Lagrangian  formulation;  we  present  its  details  below. 
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2.2.  WF  propagation 

We  implement  the  notion  of  a  WF  as  a  collection  of  points  (vertices)  forming  a  well-defined  trian¬ 
gulated  surface,  specified  by  a  mesh  (a  set  of  faces).  With  each  vertex  r  we  associate  a  ray,  which 
carries  with  itself  information  on  the  main  curvature  values  k1  =  1/ Pi  and  /c2  =  1/^2  an<^  011  a 
triad  of  vectors:  the  normal  n  to  the  WF,  and  the  directions  and  a2  of  the  main  curvatures. 
Thus,  each  vertex  on  the  WF  provides  a  local  definition  of  a  second-order  surface  (Fig.  1). 


Fig.  1:  A  local  quadratic  surface  representation  associated  with  a  ray. 

The  crucial  element  of  our  WF  implementation  is  that  the  number  of  vertices  (rays)  on  the 
WF  is  not  fixed,  but  may  change  during  its  evolution,  i.e.,  some  rays  may  be  removed  and  some 
new  rays  created.  The  essence  of  the  algorithm  is  that  the  ray-ray  distances  are  kept  in  a  certain 
range  [frmjni^maxL  controlled  by  two  parameters:  a  fixed  default  ray-ray  spacing  and  a  local 
WF  mean  curvature  radius  (defined  as  p  —  2/(1/^  |  +  |ft2|)).  Using  these  parameters,  our  evolution 
algorithm  defines,  at  each  point  of  the  WF,  minimum  and  maximum  values  of  the  ray-ray  spacing 
as 


hmin  =  7_1  min{/i0, /3  xp}  ,  (2-la) 

Vax  =  7  min{/i0,/3-1p}  ,  (2-lb) 

where  B  >  1  and  7  >  1  are  constants. 1  Eqs.  (2.1)  imply  that,  as  long  as  the  curvature  radius  p 
is  not  too  small  compared  to  h0  (p  >  /?  h0),  the  range  of  the  ray-ray  distances  is  determined  by 
the  parameter  h0;  whereas  for  smaller  curvature  radii,  the  required  values  of  ray-ray  spacing  are 
reduced  to  values  proportional  to  p.  This  prescription  helps  to  maintain  a  sufficient  resolution 
when  the  WF  shrinks  in  the  neighborhood  of  caustics.  The  value  of  the  default  resolution 
should  be  determined  on  the  basis  of  the  scale  of  the  geometrical  details  of  the  scattering  object. 

As  the  two  mesh  modification  mechanisms  that  allow  us  to  maintain  the  required  ray-ray 
spacing  we  use  the  well-known  edge  contraction  and  edge  splitting  algorithms  (Figs.  2a  and  2b);  these 
operations  enable  us  to  realize  mesh  simplification  and  mesh  “complexification”  operations.  Similar 
operations  (and  particularly  mesh  simplification)  have  been  much  studied  and  are  well  known  in 
computer  graphics.  However,  our  algorithm  differs  from  those  methods  in  one  important  aspect:  in 
creating  new  rays  (vertices)  we  are  using  information  on  the  local  second-order  surfaces  associated 
with  the  existing  rays.  This  enables  us  to  maintain  accuracy  in  the  WF  definition  of  the  second 

1  In  the  current  implementation  we  take  (3  =  4  and  7  =  2. 
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order  in  the  ray-ray  spacing.  Another  difference  between  our  mesh  simplification/ complexification 
algorithm  and  those  algorithms  commonly  used  in  computer  graphics  is  that,  at  each  contraction 
or  splitting  step,  we  sort  the  edges  according  tho  their  lengths  and  apply  contraction/splitting  to 
the  currently  shortest/longest  edge.  This  prescription  contributes  greatly  to  the  good  quality  of 
the  generated  meshes  (favorable  aspect  ratios  of  the  triangles). 


Fig.  2a:  Contraction  of  a  short  edge  to  Fig.  2b:  Splitting  of  a  long  edge  in  a  WF 
a  single  vertex  in  a  WF  mesh.  mesh. 

As  an  example,  we  show  in  Fig.  3  the  free-space  evolution  of  a  hemispherical  WF,  of  initial 
radius  R0  =  1,  discretized  with  the  resolution  h0  =  0.1.  The  WF  first  shrinks,  and  then  expands, 
with  the  evolution  step  0.3  (in  the  length  units).  Three  stages  are  shown:  step  0  (the  initial 
shrinking  WF),  with  the  radius  1.0  and  931  rays;  step  4  (just  after  the  beginning  of  the  expansion), 
with  the  radius  0.2  and  241  rays;  and  step  20,  with  the  radius  5.0  and  12,213  rays.  A  visual 
inspection  of  Fig.  3  suggests  that  a  good  quality  of  the  mesh  defining  the  WFs  is  maintained  in 
spite  of  the  wide  range  of  variation  in  the  number  of  rays. 

A  more  quantitative  assessment  of  the  algorithm  accuracy  is  presented  in  Figs.  4  and  5,  where 
deviations  of  the  WFs  from  the  ideal  spherical  shape  are  shown,  as  r.m.s.  errors.  Fig.  4  shows 
that,  in  the  later  evolution  steps,  the  errors  remain  practically  constant.  They  also  decrease  with 
the  decreasing  resolution  /i0,  although  slower  than  proportionally  to  h^.  The  latter  feature  is  due 
to  some  loss  of  accuracy  associated  with  the  WF  shrinkage  and  the  decrease  in  the  number  of  rays. 
In  a  steady  WF  expansion,  however,  the  error  is,  in  fact,  of  the  order  0(/ip)*  This  is  demonstrated 
in  Fig.  5,  where  we  plot  and  the  r.m.s.  error  in  the  surface  definition  for  a  hemispherical  WF 
expanding  in  13  steps  from  the  radius  1.0  to  the  radius  4.9.  Fig.  5  shows,  in  particular,  that  with 
decreasing  the  resolution  hQ  by  the  factor  of  2,  the  error  in  WF  definition  decreases  approximately 
by  the  factor  of  4. 
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Fig.  3:  Propagation  of  a  hemispherical  WF,  first  shrinking  (steps  0-3)  and  then  expanding  (steps 
4-20).  For  clarity,  only  halves  of  the  hemispherical  surfaces  are  shown. 


Fig.  4:  R.m.s.  error  in  the  WF  sur¬ 
face  as  a  function  of  the  evolution  step 
for  a  shrinking  and  subsequently  ex¬ 
panding  WF,  for  three  resolutions  ho  — 
0.2,0.1,0.05. 


Fig.  5:  R.m.s.  error  in  the  WF  surface 
as  a  function  of  the  evolution  step,  for 
a  steadily  expanding  WF,  for  three  res¬ 
olutions  h0  —  0.2,0.1,0.05. 


The  results  shown  above  indicate  that  we  are  able  to  change,  in  a  wide  range,  the  number  of  rays 
in  contracting  and  expanding  WFs  while  keeping  WF  definition  errors  under  control.  Such  a  precise 
interpolation  of  WF  surfaces  is  possible  only  because  in  mesh  simplification  and  complexification 
steps  we  are  using  information  on  the  local  second-order  surfaces  associated  with  rays. 


We  stress  that,  as  already  suggested  by  the  results  of  Figs3  and  4,  our  WF  discretization 
criteria  allow  for  WFs  passing  through  caustics  with  virtually  no  loss  of  accuracy. 

A  similar  example  (Fig.  6)  shows  an  ellipsoidal  WF  shrinking  and  then  expanding,  and  assum- 

ing  intermediate  shapes  of  “swallow-tail”  type. 


Fig.  6:  Evolution  of  a  WF  of  ellipsoidal  shape,  first  shrinking  and  then  expanding,  showing  creation 
of  “swallow-tail”  shapes.  (The  semitransparent  ellipsoid  in  every  evolution  step  shows  the  initial 

WF.) 

2.3.  WF  reflection  processes 

We  first  describe  our  WF  reflection  algorithm  for  smooth  scatterer  surfaces.  In  this  case  the  WF  is 
generally  a  non-smooth  surface,  i.e.,  the  tangentials  to  the  WF  may  be  discontinuous.  However, 
the  WF  remains  connected,  and  does  not  separate  into  disjoint  components. 

We  consider  reflection  of  an  incident  WF  <f>  on  an  obstacle  formed  by  a  closed  smooth  surface 
S,  as  visualized  in  Figs.  7a  and  7b.  In  this  configuration  the  initial  WF  <f>,  shown  in  Fig.  7a  as  the 
dashed  line,  approaches  the  scatterer’s  surface. 
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Fig.  7a:  A  schematic  representation  of 
one  evolution  step  of  a  WF  <j>  reflecting 
off  a  closed  impenetrable  surface  S ,  and 
forming  a  new  WF  </>'. 


Fig.  7b:  A  schematic  representation  of 
an  evolution  step  following  that  shown 
in  Fig.  7a.  SB  stands  for  the  shadow 
boundary. 


The  steps  of  the  algorithm  for  the  WF  evolution  are  then  as  follows: 

1.  We  initially  evolve  the  WF  ignoring  the  scatterer,  i.e.,  we  move  each  vertex  v  on  the  WF  <f>  by 
the  distance  As  in  the  direction  normal  to  the  WF  surface. 

2.  We  then  examine  all  vertices  v'  on  the  evolved  WF,  and  classify  them  into  “inner”  and  “outer” 
ones:  inner  vertices  v'  are  such  that  the  ray  segment  (vv1)  intersects  the  surface,  and  outer 
vertices  are  the  remaining  ones.  In  Fig.  7a  we  marked  inner  vertices  with  red  dots,  and  outer 
vertices  with  black  dots. 

3.  Next,  we  examine  all  edges  of  the  evolved  WF.  If  an  edge  connects  an  inner  vertex  and  an 
outer  vertex,  such  as  the  vertices  a'  and  b'  in  Fig.  7a,  we  find,  by  performing  a  binary  search, 
the  point  where  it  intersects  the  surface  S.  At  this  point  of  the  edge  we  create  two  new  vertices, 
an  “inner”  and  an  “outer”  one.  Such  vertices  in  Fig.  7a  are  a[  and  a!0.  They  share  the  same 
spatial  location,  but  will  follow  different  evolution  paths. 

4.  We  examine  in  turn  all  faces  (triangles)  of  the  evolved  WF  mesh,  and  process  “divided” faces, 
i.e.,  those  faces  whose  some  vertices  are  inner  and  others  are  outer.  An  example  of  such  a  face 
is  shown  in  Fig.  8.  We  partition  such  a  face  into  a  triangle  (an  outer  triangle  in  Fig.  8)  and  a 
quadrilateral,  and  split  further  the  quadrilateral  into  two  triangles. 

5.  Finally,  we  reprocess  all  the  inner  vertices  (the  original  evolved  vertices  and  the  new  ones).  We 
evolve  them  back  to  the  original  WF  <j>,  and  evolve  them  again  forward  in  the  presence  of  the 
scatterer ,  i.e.,  taking  reflections  into  account.  In  this  way  the  inner  vertices  form  a  part  of  the 
reflected  WF  <f>'  (the  red  dash-dotted  line  in  Fig.  7a),  while  the  previously  found  outer  vertices 
form  the  forward-propagating  part  of  the  WF  <j>'  (the  black  dash-dotted  line  in  Fig.  7). 
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Fig.  8:  Partition  of  a  WF  face  having  two  inner  vertices  {a!  and  c')  and  one  outer  vertex  ( b '). 

We  stress  an  important  element  of  our  algorithm:  The  mesh  faces  such  as  shown  in  Fig.  8  are 
divided  into  parts  belonging  to  the  forward-propagating  and  reflected  segments  of  the  WF.  The 
algorithm  is  able  to  locate  the  division  points  (new  vertices  a[,  a'Q,  7-,  and  7'  in  Fig.  8)  with  an 
arbitrary  accuracy ,  unrelated  to  the  ray-ray  spacing  h.  Therefore,  the  location  rays  on  the  boundary 
between  the  forward-propagating  and  reflected  parts  of  the  WF  can  be  considered  exact,  the  only 
approximation  made  is  the  considered  boundary  consists  of  straight  line  segments,  rather  than 
curved  lines.  For  typical  segment  length  of  order  h0,  the  resulting  error  is  0(hl/R),  where  R  is 
the  curvature  radius  of  the  surface  S. 

The  algorithm  described  above  can  be  generalized  to  scatterers  formed  by  non-smooth  surfaces, 
both  closed  and  open.  The  generalization  is,  however,  somewhat  more  involved,  for  two  main 

reasons: 

(i)  the  WF  usually  separates  into  several  disconnected  pieces;  and 

(ii)  in  general,  it  is  necessary  to  consider  more  complex  cases  of  the  scatterer  surface  intersecting 
the  WF  mesh;  in  particular,  triangles  forming  the  WF  mesh  may  be  intersected  in  ways  more 
complicated  than  that  shown  in  Fig.  7  (in  the  analysis  of  smooth  surfaces  we  assumed  that 
the  surface’s  radii  of  curvature  R  were  much  larger  than  the  ray-ray  spacing  h). 

If  we  attempt  to  reproduce  GO  shadow  and  reflection  boundaries  with  an  arbitrary  accuracy,  the 
second  of  these  problems  causes  a  significant  complication  of  the  algorithm.  Therefore,  at  the 
present  stage  of  the  algorithm  development,  we  concentrated  on  the  first  problem,  that  is,  splitting 
of  the  WF  due  to  reflection  on  “sharp”  edges  of  the  scatterer. 

We  give  below  a  brief  description  of  the  WF  reflection  algorithm  for  a  practically  important 
case  of  a  scatterer  with  a  facetized  (typically  triangulated)  surface  representation.  We  consider  a 
scatterer  whose  part  is  represented  schematically  in  Fig.  9.  The  initial  incident  WF,  denoted  by 
<j>,  approaches  a  wedge-type  element  of  the  scatterer  surface  S.  The  edge  of  the  wedge  shown  in 
the  figure  is  assumed  to  be  “sharp”,  i.e.,  the  angle  between  the  normals  to  the  faces  exceeds  some 
predefined,  critical  angle  60. 

As  before,  we  first  evolve  the  rays  associated  with  the  WF  as  if  the  scatterer  were  absent; 
subsequently,  we  examine  intersections  of  rays  with  the  scatterer  surface,  and  then  create  some 
additional  rays,  and  reflect  ray  segments  penetrating  the  scatterer  surface. 
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Fig.  9a:  A  schematic  representation  of  Fig.  9b:  The  same  WF  configuration  as 
one  evolution  step  of  a  WF  (j)  reflecting  in  Fig.  9a,  with  the  rays  and  vertices 
off  a  non-smooth  impenetrable  surface  removed  for  clarity. 

S',  and  forming  a  new  WF  ft. 

More  specifically,  the  main  steps  of  the  algorithm  are  as  follows: 

1.  We  evolve  the  WF  ignoring  the  scatterer,  i.e.,  we  move  each  vertex  on  the  WF  <f>  by  the  distance 
As  in  the  direction  normal  to  the  WF  surface;  thus  the  vertices  a,  6, . . .  become  a',  b\ . . .  . 

2.  We  locate  sharp  edges  on  the  surface  S  -  in  our  example  the  tip  of  the  wedge.  This  task  is 
carried  out  as  follows: 

2.1.  For  each  ray  segment  intersecting  the  surface  S  we  find  the  normal  to  the  surface  at  the 
intersection  point. 

2.2.  For  each  edge  on  the  WF  mesh  we  compare  normals  at  the  two  ends.  For  example,  for  the 
original  WF  mesh  edge  (a  b)  we  compare  the  normals  at  points  where  the  ray  segments 
(a a')  and  ( bb ')  intersect  S. 

2.3.  If  the  angle  between  the  normals  at  the  end  vertices  of  an  edge  exceeds  the  critical  angle, 
we  examine  the  edge  in  more  detail:  we  subdivide  it  in  a  binary  fashion  until  we  find  two 
ray  segments  intersecting  the  surface  S  very  near  the  edge  (this  can  be  done  with  any 
predefined  accuracy).  In  Fig.  9  these  rays  are  (51  <5^),  and  (ft  ft). 

3.  As  in  the  previous  algorithm,  we  identify  then  pairs  of  rays  with  end  slightly  outside  and 
slightly  inside  the  scatterer  (again  with  some  predefined  accuracy).  Such  ray  pairs  in  Fig.  9 
are  {(o^  ft) ,  (a2  ft)}  and  {(ft  ft) ,  (ft  ft)}* 

4.  We  split  (partition)  faces  whose  edges  contain  the  additional  vertices  associated  with  ray 
segments  created  in  the  previous  step. 

5.  Finally,  as  before,  we  reflect  the  rays  which  have  penetrated  the  scatterer  surface  S . 

In  the  example  illustrated  in  Figs.  9  the  resulting  reflected  WF  consists  of  several  segments 
(Fig.  9b):  the  parts  ft  and  ft  of  the  incident  WF  which  have  not  yet  reached  the  scatterer’s 
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surface,  and  the  parts  (f>'3  and  4>'i  which  have  been  reflected  from  the  two  facets  forming  the  wedge. 
Reflection  on  a  surface  with  an  edge  creates  thus  a  “hole”  in  the  WF,  unlike  reflection  on  a  smooth 
surface  (Figs.  7a  and  7b),  where  the  WF  remains  a  connected  surface.  (We  note  that,  as  we  will 
see,  such  “holes”  are  filled  in  by  WF  due  to  diffraction.) 

As  an  example  of  WF  evolution  in  the  presence  of  smooth  scatterers,  we  consider  a  process 
of  multiple  reflection  on  a  system  of  two  spheres,  S1  and  S2  (Fig.  10).  The  radii  of  the  spheres 
are  equal  1,  and  the  wavefront  resolution  is  h0  =  0.1.  The  initial  planar  wavefront  travels  upward 
at  at  angle  of  approximately  45°.  The  consecutive  evolution  steps  show  reflections  from  the  first, 
second,  and  again  from  the  first  sphere. 


step  4  «»ep5  «*eP6 


Fig.  10:  Reflection  of  a  planar  WF  on  a  system  of  two  spheres. 

In  Fig.  11a  we  show  an  evolution  stage  for  the  same  system  of  two  spheres  arrived  at  after 
several  reflections.  (In  this  example,  the  initial  planar  wavefront  travels  downward  in  the  vertical 
direction.) 


12 


Fig.  11a:  Reflection  of  a  planar  WF  on 
a  system  of  two  spheres,  for  the  resolu¬ 
tion  h0  =  0.1,  and  the  number  of  rays 
increasing  from  1,681  to  22,962.  Var¬ 
ious  colors  show  the  parts  of  the  WF 
W0, ...  ,W3  due  to  0,  1,  2,  and  3  reflec¬ 
tions. 


Fig.  lib:  The  same  final  WF  as  shown 
in  Fig.  11a,  but  obtained  using  the  con¬ 
ventional  ray-tracing  algorithm  with  a 
much  higher  resolution  h0  =  0.02  and  a 
fixed  number  of  rays  40,401. 


Parts  of  the  wavefront  due  to  0,  1,  2,  and  3  reflections  are  shown  in  various  colors.  All  these 
wavefront  parts  are  equally  well  defined,  due  to  the  dynamical  ray  number  adjustment.  In  the 
considered  process  of  wavefront  scattering  and  expansion  the  number  of  rays  increases  from  N  = 
1,681  to  N  =  22,962. 

For  comparison,  we  show  in  Fig.  lib  the  analogous  final  wavefront  constructed  using  the 
conventional  ray-tracing  algorithm  with  a  much  higher  resolution,  h0  —  0.02,  and  the  corresponding 
fixed  number  of  rays  NRT  =  40,401.  This  Figure  clearly  shows  the  rapid  deterioration  of  the 
wavefront  definition  with  the  increasing  reflection  order  and  the  associated  wavefront  expansion. 
In  particular,  the  last  wavefront  IV3,  due  to  triple  reflection,  is  barely  existent.  The  comparison 
indicates  that,  due  to  the  dynamic  ray  number  adjustment,  our  wavefront  evolution  algorithm 
provides  a  much  superior  resolution  at  all  reflection  orders  with  an  overall  much  smaller  number  of 
rays. 


2.4.  WF  edge  diffraction  processes 

We  describe  now  our  WF  algorithm  for  edge  diffraction  processes.  While  the  algorithm  is  already 
complete,  we  are  currently  in  the  process  of  its  numerical  implementation,  and  thus  the  results 
presented  here  are  preliminary. 

The  geometrical  construction  of  edge  diffracted  WFs  is  based  on  the  original  Geometrical 
Theory  of  Diffraction  (GTD)  [12, 13],  which  in  turn  follows  from  the  generalized  Fermat’s  principle, 
and  ultimately  from  the  stationary-point  asymptotics  of  the  high-frequency  scattering  processes. 
The  associated  field  amplitudes,  however,  are  computed  using  the  Uniform  Geometrical  Theory  of 
Diffraction  (UTD)  [14,15]. 

The  main  steps  of  the  algorithm  for  construction  of  edge-diffracted  WFs  are  as  follows: 

1.  Identify  scatterer  edges  constituting  source  of  diffraction  (we  refer  to  these  as  “diffraction 
edges”). 


2.  Generate  diffracted  rays  emanating  from  the  diffraction  edges. 

3.  Assign  WF  curvature  parameters  and  field  values  to  the  created  diffracted  rays. 

4.  Construct  triangulation  mesh  (connectivity)  for  the  diffracted  rays. 

Although  the  above  steps  may  seem  straightforward,  their  implementation  is  rather  complex.  We 

discuss  them  now  in  turn. 

Step  1  requires  implementing  an  algorithm  with  elements  of  the  previously  discussed  proce¬ 
dures  for  (1)  detecting  “sharp  edges”  (in  the  context  of  WF  splitting  associated  with  reflection), 
and  (2)  detecting  the  shadow  boundary  (described  in  the  context  of  the  smooth-surface  reflection 
algorithm).  If  the  edge-diffraction  algorithm  is  applied  to  a  facetized  scatterer  surface,  we  assume 
that  the  “sharp  edges”  are  sources  of  diffraction,  while  “smooth  edges”  have  to  the  treated  by  a 

separate  algorithm  for  smooth-surface  diffraction. 

An  important  element  of  the  algorithm  is  that  the  identified  “diffraction  edges”  are  associated 
with  and  accessed  through  the  WF  mesh  edges.  As  indicated  in  Fig.  12,  new  rays  (a1?a2)...  ) 
hitting  the  sharp  edge  Eon  the  scatterer  emerge  from  the  edges  of  the  WF  mesh  {ex,e2,  ■  ■  ■  )• 
Indeed,  these  rays  are  created  by  interpolating  rays  at  end  vertices  of  the  WF  edges. 


Fig.  12:  A  schematic  view  of  the  procedure  of  generating  diffracted  rays  emerging  from  the  points 
of  a  diffraction  edge  E  hit  by  newly  created  rays  ax,a2,.... 

Step  2  is  performed  in  a  loop  through  the  edges  e  of  the  WF .  From  some  of  these  edges  there 
emerges  a  pair  of  new  rays  a,  terminating  at  a  diffraction  edge.  At  the  termination  point  we  create 
a  set  of  new  rays  (of  length  equal  to  the  evolution  step  decreased  by  the  length  of  the  ray  segments 
a),  on  the  surface  of  a  cone  (the  “Keller  cone”),  whose  axis  is  the  diffraction  edge  (Fig.  12).  The 
number  of  diffracted  rays  is  fixed,  and  is  determined  on  the  basis  of  the  required  WF  resolution 
(this  element  of  the  algorithm  is  not  critical,  since  the  number  of  rays  will  adjust  itself  during  the 
subsequent  evolution  of  the  diffracted  WF). 

Finally,  step  3  is  executed  as  a  loop  over  the  faces  F  of  the  WF.  We  consider  here  such  faces 
F  from  two  edges  of  which  emerge  new  rays  a  (Fig.  13a.)  Through  these  rays  we  access  two 
consecutive  sets  of  diffracted  rays,  and  create  edges  connecting  the  corresponding  rays  in  the  two 
sets,  forming  a  surface  consisting  of  quadrilateral  facets.  By  triangulating  that  surface  (Fig.  13b), 
we  arrive  at  a  meshed  diffracted  WF . 
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Fig.  13a:  Building  connectivity  for  the  Fig.  13b:  TViangulation  of  the  diffracted 
diffracted  rays.  WF. 

The  overall  picture  of  the  newly  created  diffracted  WF  may  be  as  shown  in  Fig.  14.  In  this 
Figure  a  plane  WF  <j>  is  approaching  a  smooth  edge  of  a  screen  S,  moving  in  a  direction  nearly 
tangential  to  the  screen  plane.  As  the  incident  WF  advances  to  the  position  </>',  it  generates  a 
reflected  WF  (not  shown),  and  a  diffracted  WF,  according  to  the  algorithm  described  above. 


Fig.  14a:  A  schematic  view  of  a  new  Fig.  14b:  A  perspective  view  of  the 

diffracted  WF  generated  along  a  smooth  diffracted  WF  shown  in  Fig.  14a. 

edge  of  a  screen  <5  by  a  plane  WF  inci¬ 
dent  nearly  tangentially  to  the  screen. 

We  note  that  the  diffracted  WF  algorithm,  as  described  above,  may  not  reproduce  correctly 
the  end  (cone-shaped)  parts  of  the  new  WF;  this  effect  is  indicated  in  Fig.  15,  which  shows  the 
actual  process  of  diffracted  WF  creation,  as  modeled  by  our  code.  The  reason  for  this  deficiency  is 
that  we  use  only  WF  edges  to  detect  WF  intersections  with  diffraction  edges.  In  other  words,  the 
pairs  of  new  rays  a,  terminating  at  the  diffraction  edge  E,  may  emerge  only  from  the  WF  edges 
g  'j'jiig  feature  of  our  procedure  affects  also  the  GO  reflection  algorithm,  where  some  parts  of  the 
reflected  WF  may  be  missing.  In  both  cases  the  defects  in  the  WFs  are  only  due  to  the  present 
implementation  of  the  algorithm,  and  will  be  corrected  as  the  algorithm  is  further  refined. 


15 


diffracted  WF  segment  diffracted  WF  segment 
created  in  step  2  created  in  step  1 


WF  defect 


Fig  15:  A  defect  of  the  created  WF  -  Fig.  16:  A  gap  between  the  diffracted 

a  missing  end  of  the  cone-shaped  WF  WF  segments  created  in  two  consecutive 

part,  pointed  to  by  the  arrow.  evolution  steps. 

Another,  related,  defect  in  the  present  algorithm  for  diffracted  WF  creation  manifests  itself 
in  gaps  between  the  segments  of  the  diffracted  WF  generated  in  consecutive  evolution  steps.  A 
situation  of  this  type  is  shown  for  diffraction  on  a  disc  in  Fig.  16.  These  deficiencies  can  be  corrected 
by  extending  the  “memory”  of  the  algorithm,  i.e.,  increasing  the  amount  of  information  preserved 
from  one  to  the  next  evolution  step. 

The  present  version  of  our  code  includes  only  edge  diffraction  as  an  asymptotic  expansion 
scattering  mechanism  of  order  1/Vk.  Another  diffraction  mechanism  not  implemented  in  the 
current  version  of  our  code  is  diffraction  from  corners,  inevitably  arising  for  diffraction  edges 
consisting  of  straight  line  segments.  As  an  example,  we  show  in  Fig.  17  several  steps  of  creating 
diffracted  WFs  due  to  a  plane  WF  incident  on  a  rectangular  plate.  The  WF  diffracted  on  the 
individual  edges  of  the  plate  form  either  cylinders  of  cut-off  cones.  These  WFs  should  be  connected 
by  the  four  spherical  WFs  due  to  diffraction  on  the  four  corners  of  the  plate. 
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scatterer  surface  (plate) 


Fig.  17:  Three  evolution  steps  in  creating  diffracted  WFs  due  to  a  plane  WF  incident  on  a  rect¬ 
angular  plate.  For  clarity,  the  primary  and  reflected  WFs  are  not  shown  in  evolution  steps  other 
than  the  initial  configuration  ( “step  0” ) . 


3.  Field  computation 

3.1.  Behavior  of  edge  diffracted  fields 

The  final  outcome  of  the  WF  algorithm  should  be  in  the  form  of  the  observable  quantities  - 
typically  the  scattered  fields,  and  the  related  scattering  cross-section.  Various  approaches  to  this 
goal  are  possible,  and  we  consider  two  of  them: 

1.  The  scattered  fields  can  be  determined  directly  from  the  fields  associated  with  the  sufficiently 
evolved  WFs,  i.e.  WFs  leaving  the  scattering  region. 

Alternatively, 

2.  The  scattered  fields  can  be  computed  as  Kirchhoff  integrals  of  the  currents  induced  on  the 
scatterer  surface  (more  generally,  effective  currents  on  an  appropriately  defined  aperture). 

The  two  approaches  are  equivalent,  provided  all  terms  to  the  given  order  in  1/k  in  the  asymptotic 
expansion  are  consistently  taken  into  account,  both  in  computation  of  the  currents  and  computation 
of  the  field.  In  both  cases  we  utilize  fields  associated  with  rays  constituting  WFs.  In  particular, 
the  second  approach  requires  the  electric  induced  current 

J  =  -2n  x  AH  ,  (3-1) 
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where  AH  is  the  discontinuity  of  the  magnetic  field  across  the  scatterer  surface,  and  n  is  the 
normal  to  the  surface.  Because  of  this  fact,  we  find  it  more  convenient  to  use  the  magnetic  rather 

than  electric  field  as  the  basic  field  associated  with  rays. 

Thus,  to  each  point  of  a  WF  we  ascribe  a  magnetic  field  value  H;  a  field  value  is  also  associated 
with  every  point  on  a  ray  connecting  two  consecutive  WFs.  For  a  WF  propagating  in  free  space, 
knowing  the  field  value  H(s0)  at  a  certain  evolution  parameter  s0,  we  can  obtain  (according  to  the 
laws  of  GO)  the  field  value  H(s)  at  another  value  s  of  the  evolution  parameter  as 


H  (s) 


P1P2 


( Pi  +  s  ~  »o)  (^2  "b  s  so) 


H(s0) 


(3.2) 


where  the  curvature  radii  refer  to  the  WF  at  the  reference  evolution  parameter  s0.  Here  the  fields 
H  are  defined  so  as  to  exclude  the  common  phase  factor  exp(ifcs)  associated  with  the  evolution 

parameter  s. 

When  creating  new  rays  due  to  mesh  complexification  or  simplification,  we  associate  with  them 
field  values  obtained  by  interpolation  of  the  field  values  of  the  neighboring  rays.  In  the  present 
implementation,  when  creating  new  ray  in  edge  contraction  (Fig.  2a)  or  edge  splitting  (Fig.  2b)  we 
use  simple  linear  interpolation  of  field  values  between  the  end  vertices  of  the  considered  edge. 

When  creating  diffracted  rays,  we  ascribe  to  them  fields  calculated  according  to  UTD  [14, 15], 
or  rather  to  its  modification  [1],  which  we  present  below.  Our  modification  of  UTD  was  originally 
motivated  by  the  need  of  calculating  induced  currents  near  the  diffraction  sources  (edges);  we  thus 
refer  to  it  as  nUTD,  or  the  “ near-field  UTD”.  We  also  discuss  it  mostly  in  the  context  of  currents, 
although  similar  considerations  apply  to  the  general  problem  of  field  computation. 

It  is  important  to  stress  that,  in  the  “transition  regions”  (i.e.  near  the  diffraction  source  and 
near  the  optical  shadow  and  reflection  boundaries),  the  fields  provided  by  UTD  or  nUTD  are 
not  ray-optical ,  i.e.,  they  do  not  evolve  according  to  Eq.  (3.2),  but  are  given  by  more  complex 
expressions  involving  appropriate  transition  functions.  Only  at  large  distances  from  the  diffraction 
source,  and  outside  the  transition  regions,  do  they  reduce  to  ray-optical  fields  of  GTD. 

During  Phase  I,  we  addressed  the  problem  of  modifying  fields  near  the  diffraction  sources 
only  for  diffraction  on  a  flat  polygonal  screen  with  straight  edges,  illuminated  with  a  plane  wave. 
Near  the  edges  of  a  polygonal  screen,  the  fields  are  given,  asymptotically,  by  the  solution  of  the 
corresponding  canonical  problem  -  in  this  case,  a  half-plane.  Hence,  we  first  examine  the  behavior 
of  the  fields  and  currents  obtained  form  the  exact  half-plane  solution.  Then,  assuming  the  exact 
half  plane  solution  for  the  field  diffracted  by  the  edges  of  the  screen,  we  compute  the  current 
induced  on  the  screen. 

Behavior  of  the  fields  and  currents  near  diffraction  edges 

Consider  a  half-plane  defined  by  0  <  x  <  00,  y  =  0,  —00  <  z  <  00  (Fig.  18).  We  assume  an 
incident  plane  wave  with  the  wave  vector 


k;  =  —  k  (sin  /?0  cos  </>0  x  +  sin/30  sin  <p0  y  +  cos  /30  z) 


(3.3) 


(incidence  normal  to  the  edge  would  correspond  to  Pq  —  tt/2). 
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Incident  wave  direction 


half-plane 

Fig.  18:  The  half-plane  geometry  and  the  field  configuration  for  the  “H-polarization  problem. 

The  exact  solution  for  the  canonical  problem  is  well  known  (see,  e.g.,  Ref.  16]).  The  Maxwell 
equations  separate  in  this  case  into  independent  Helmholtz  equations  for  the  z-components  of  the 
electric  and  magnetic  fields,  with  appropriate  boundary  conditions. 

We  first  consider  the  “H-polarization”  problem  for  the  magnetic  field  Hz  (the  component  along 
the  edge  of  the  half-plane).  The  solution  contains  the  incident,  reflected,  and  diffracted  fields,  the 
latter  being 

Hz(p,  <t>)  =  — eiir/4  eikp/sin0o  {sec  F  (2*P  sin  A)  cos2  ^2^) 

+  sec  -  +^—  F  ^2  kp  sin  j30  cos2  j  H'z  , 

(3.4) 


where  the  function  F,  related  to  the  Fresnel  integral, 

OO 

F(x)  =  -2iv/xe_il  J  dfe5*2  , 

VS 


is  normalized  so  that 


and 


F(x)  oc  1  -  —  +  •  •  •  for  x  ->•  oo  , 

JLtX 


F(x)  =  [V^  -2xe~i*'4  +  ---]e-i  <*+,r/4)  for  x  ->  0  , 
and  satisfies  the  relation  . 

The  discontinuity  of  the  magnetic  field  across  the  screen  is  then 

AHdz(p)  =  H*(p,0)-H«(p,2ir) 

_  _ei  w/4  ei  fcp/  sin  0O  2  \/s^  A)  sec  ^2.  F  (2  kp  sin  /30  cos2  , 
y/2 nkp  2  \  2  / 


(3.5) 


(3.6) 

(3.7) 

(3.8) 


(3.9) 
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(3.10) 


and,  according  to  Eq.(3.1)  (with  n  =  y),  the  electric  current  due  to  diffraction  is 

4(p)  =  -2 A HAZ{P)  • 

As  follows  from  Eq.(3.6),  for  large  distances  from  the  edge  (i.e.,  kp  cos2 (4/2)  »  1),  the 
current  behaves  like  const/ y/kp,  which  is  a  reflection  of  the  cylindrical- wave  nature  of  the  diffracted 
fields.  Near  the  edge,  on  the  other  hand,  Eq.(3.7)  implies  that  Jd(p)  tends  to  a  constant,  which, 
actually,  cancels  the  PO  contribution  from  the  current  terms  j£(p)  and  Jx(p)  due  to  the  incident 
and  reflected  waves.  A  more  detailed  analysis  (using  the  expansion  (3.7)  for  x  ->•  0)  shows  that 
the  total  current  exhibits  the  correct  edge  behavior,  Jx(p)  ~  y/P  for  p  -4  0. 

Similar  analysis  applies  to  the  “E-polarization”  problem,  in  which  the  ^-component  of  the 

diffracted  electric  is 

E*M)  =  -e,7r/4eifcp/sin^°  7Q§=p  |sec^2F  (ikp  sinfo  cos2 

-  sec  F  ^2  kp  sin  4  cos2  t+jtsi'j  |  Elz  . 

(3.11) 

Here,  after  calculating  the  magnetic  field  and  its  discontinuity,  we  find  that  the  total  current  (now 
directed  along  the  z-axis)  is  Jz(p)  ~  1/y/p  for  r  ->  0,  again  in  agreement  with  the  expected  edge 
behavior. 

Behavior  of  conventional  TJTD  fields  and  currents  near  diffraction  edges 
Having  established  the  behavior  of  the  exact  solution  for  the  canonical  problem,  we  can  now 
compare  them  with  the  fields  predicted  by  UTD  [14]  (converted  to  our  conventions,  and  specialized 
to  the  case  of  a  screen  with  a  straight  edge  illuminated  by  a  plane  wave).  As  we  will  see,  these 
expressions  will  require  some  modification  of  the  conventional  UTD  formulation  in  order  to  correctly 
reproduce  fields  and  currents  in  the  vicinity  of  the  diffraction  edge. 

To  specify  the  UTD  diffraction  coefficients,  we  first  have  to  define  the  relevant  geometrical 
quantities.  Thus,  in  the  reference  system  of  Fig.  18,  the  edge  direction  is 

e  =  z  (3-12) 

and  the  directions  of  the  incident  and  diffracted  rays  are  given  by 

s1  =  -  sin 4  cos 4  x  -  sin 4  sin 4  y  -  cos 4  z  ,  sd  =  sin/3  cos </>  x  +  sin/3  sin </>  y  -  cos/3  z  , 

(3.13) 

where  the  diffracted  ray  is  required  to  lie  on  the  diffraction  cone,  i.e.,  e  •  sd  =  e  •  s1,  or  /3  =  4- 
These  two  vectors  and  the  vector  e  define  two  planes:  the  edge-fixed  plane  of  incidence  (e,  s  ),  and 
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the  edge-fixed  plane  of  diffraction  (e.  sd ) .  The  unit  directional  vectors  associated  with  the  angles 


(j> o,  <f>,  P0,  and  p  are  then 

<j>0  =  -  sin  <t> o  x  +  cos  4>q  y  ,  (3.14a) 

0  =  —  sin  (f>  x  +  cos  <f>  y  ,  (3.14b) 

p0  =  cosfio  cos  d>0  x  +  cospQ  sin  <j>0  y  -  sin/30  z  ,  (3.14c) 

£  =  cos  /30  cos  <j)  x  +  cos  /30  sin  <f>  y  +  sin/30  z  .  (3.14d) 


It  is  now  convenient  to  express  the  incident  and  diffracted  fields  in  terms  of  the  triads  of  unit 
vectors  (s '\<j>0iPo)  an<l  (sd,<£,/3)  associated  with  the  incident  and  diffracted  rays.  The  diffracted 
field  Ed  at  the  observation  point  r  is  then 

Ed(uTD)(r)  _  *  e'**"  W.  B‘(r0) 

vsd  1 

+-^h(A)i ‘An  ^  s<i)  '  ®*(ro)}  i  (3.15) 

in  terms  of  the  incident  field  E'  at  the  corresponding  diffraction  point  Tq  (located  on  the  edge), 
the  path  length  of  the  diffracted  ray  sd  =  |r  -  r0|,  and  the  “soft”  and  “hard”  UTD  diffraction 
coefficients  Ds  and  Dh.  The  latter  are  given  by 

B.(  A,  *,,**,*•)  =  {sec  ±^LF(2ks‘  +ffi,  cos’  tjl) 

=F sec  p  ^2 k sd  sin 2P0  cos2  t+jtl ^  ^  , 

(3.16) 


where  the  function  F  is,  as  before,  defined  by  Eq.(3.5). 

Since  p=  sA  sin  P0,  we  can  see  that  the  arguments  of  the  functions  F  in  Eq.(3.16)  are  exactly 
the  same  as  in  Eqs.  (3.4)  and  (3.11).  Verification  of  Eq.(3.11)  for  the  E-polarization  problem  is 
now  straightforward.  By  taking  the  ^-component  of  the  diffracted  electric  field  (3.15)  we  find  that 
only  the  soft  diffraction  coefficient  contributes,  because  z  •  0  =  0.  Since  the  E-polarization  implies 
incident  unit  field  E1  =  —fio,  we  have  Pq  •  E!  =  —1,  and  substituting  Ds  of  Eq.(3.16)  into  Eq.(3.15), 
we  immediately  reproduce  Eq.(3.11),  i.e., 

(UTD)(r)  _  Eg{r)  .  (3.17) 

In  order  to  compare  the  UTD  expression  with  the  exact  solution  (3.4)  for  the  H-polarization, 
we  first  use  the  exact  relation  (following  from  the  Maxwell  equations) 

Hd  (UTD)  =  _  1 V  X  Ed  (UTD)  ,  (3.18) 

k 

expressing  the  magnetic  field  in  terms  of  the  diffracted  electric  field  (3.16). 
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To  simplify  the  result,  we  specialize  now  to  normal  incidence,  fi0  =  7r/2.  111  this  case  we  have 


h?  <utd)  (p, « = 4 ».  <OTD  w> 


=  eiW4  ei  kp 


|  sec 


+  sec 


2  y/2%kp 

<l>-  <fto 

2 

<l>  +  <f>  o 


1  +  cos (<f>  -  (f> o) 


1  -  F  (  2  fcp  cos: 


1  +  cos((/)  +  4>q) 


1-f(\ 


2  kp  cos 


2 

2 

2  <^  +  00 


)]] 


Hi . 


(3.19) 


Eq.(3.6)  implies  that  Eq.(3.19)  approaches  the  exact  result  (3.4)  for  large  kp.  To  assess  the  signif¬ 
icance  of  the  differences  for  small  values  of  kp  for  the  computation  of  the  current  (Eq.(3.10)),  we 
evaluate  the  discontinuity  of  the  magnetic  field  across  the  screen, 


A H^VTD\p)  =  i^d(UTD)(p,0)  -  ifyd(UTD)(p,27r) 

2 


=  ei7r/4 


pifcp  *  sec  — 

■\Z2nkp  2 


jl+cos<£0  1  -  F  ^2 kp  cos2  |  Hz  • 


(3.20) 


This  results  agrees  with  the  exact  Eq.(3.9)  (for  f30  =  tt/2)  only  for  large  arguments  of  the  transition 
function  F,  but  not  for  small  arguments.  For  p  ^  0  the  UTD  result  is  actually  singular  (~  1/^/p), 
while  the  exact  expression  (3.9)  tends  to  a  constant. 

We  have  thus  found  that  if  we  use  the  UTD  formula  (3.15)  to  calculate  the  electric  field,  and 
then  evaluate  the  magnetic  field  using  the  exact  relation  following  from  Maxwell’s  equations,  the 
result  is  incorrect.  On  the  other  hand,  we  obtain  a  correct  result  if  we  express  the  magnetic  field 
in  terms  of  the  electric  field  using  the  approximate  relation 

jjd  (UTD)  ^  §d  x  Ed  (UTD)  ^  (3.21) 

valid  only  asymptotically,  and  assuming  the  ray-optical  nature  of  the  electric  field.  Eq.(3.21) 
implies  now 

tfd  (UTD)(r)  ^  z  .  (§d  x  Ed  (UTD)(r))  =  — eiks*  z  -pDh(0o,  <f>0 ,  </>,  k,  sd)  4>0  •  E^o)  ,  (3.22) 

z  vs“ 

which  coincides  exactly  with  Eq.(3.4). 

It  would,  therefore,  seem  that  UTD  does  reproduce  the  correct  fields  near  the  diffraction  edge, 
provided  the  magnetic  field  is  calculated  according  to  Eq.(3.21),  and  not  Eq.(3.18).  This,  however, 
is  not  an  acceptable  solution  for  computing  all  components  of  the  magnetic  field.  In  particular,  if 
we  consider  again  the  E-polarization  case,  taking  for  simplicity  normal  incidence  (^30  =  7r/2),  the 
electric  field  is  directed  along  the  edge  (the  z-axis)  and  is  given  by  Eqs.  (3.11)  and  (3.17).  In  this 
case,  if  we  evaluate  the  magnetic  field  on  the  surface  of  the  screen,  we  have  sd  =  x,  and  Eq.(3.21) 
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yields  the  magnetic  field  in  the  direction  normal  to  the  screen,  sd  x  Ed(UTD)  ~  y.  Therefore, 
the  resulting  electric  current  (Eq.(3.1))  is  exactly  zero ,  leading  to  the  incorrect  conclusion  that 
soft  diffraction  generates  no  scattered  field  at  all;  while  in  reality  the  current  is  nonzero,  because 
a  tangential  magnetic  field  results  from  the  normal  derivative  of  the  tangential  electric  field,  and, 
according  to  Eq.(3.18),  we  have  Hx  ^  ^  —  —i/k  dy  Ez 

Behavior  of  modified  UTD  (nUTD)  fields  and  currents  near  diffraction  edges 
We  found  above  that  the  conventional  UTD  does  not,  in  general,  reproduce  correct  near  fields  and 

thus  surface  currents  near  the  diffraction  edges. 

In  this  situation  we  suggest  a  modification  of  the  UTD  procedure  to  allow  for  a  correct  de¬ 
scription  of  all  field  components,  and  a  reliable  computation  of  surface  currents.  Here  we  sketch 
the  modified  formulation  only  for  a  straight  edge  of  a  flat  screen;  generalization  to  other  cases  is 
more  complex  but  feasible. 

Our  starting  point  are  the  exact  half-plane  solutions  (3.11)  and  (3.4)  for  the  E-  and  H- 
polarization.  These  solutions  can  be,  without  any  approximation,  represented  as 


£d  =  -4=eifcsDs4  ,  (3.23a) 

V5 

H*  =  --^eiksDhHl  ,  (3.23b) 

Vs 


where  s  =  sd. 

We  use  Eq.(3.23b)  for  the  z-component  of  the  magnetic  field,  and  express  the  remaining 
components  of  Hd  in  terms  of  tfd  and  £d,  by  applying  the  well  known  formula  for  the  transverse 
components  of  the  magnetic  field  in  terms  of  the  z-components,  valid  for  solution  of  Maxwell 
equation  having  an  exponential  z-dependence,  in  our  case  exp(— i  kz  cos/30).  We  also  express  the 
incident  fields  E\  and  Hlz  in  terms  of  the  incident  magnetic  field  components  along  the  triad  of 
vectors  (s\<p0,j30).  Assuming  transversality  of  the  incident  fields,  we  have  H1  =  <j)0 11^  +  j30Hp, 
and,  from  Maxwell’s  equations, 

E\  =  sin  /30  H]p  ,  H'z  =  -  sin  fi0H'p  .  (3.24) 


Finally,  we  substitute  in  the  resulting  expression  the  relations  Eqs.  (3.23).  In  this  way  we  find  the 
total  diffracted  field  as 


Hd  =  — [*  *  VT£d-cos/?0VTtfd]+zifd 
k  sm  /30 

= - V—  {sin/30z  x  VTDsHl4)  +  cos£0  sin/30  VxDhJ?^}  ~fzelkS 

k  sin  Pq  v 

+  z  sin^0 £>h  Hp  -^=  eifes  . 


(3.25) 


We  use  now 


VT  =  p  sin  P0ds  +  4> 


s  sin  Pq 


d, 


<P  ’ 


(3.26) 
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where 


(3.27) 


p  =  cos  <f>  x  +  sin  <j)  y  . 

To  the  leading  order,  the  gradients  in  Eq.(3.25)  act  only  on  the  factor  exp(i  ks), 

VTeifcs  =  pik  sin/30elfcs  ;  (3.28) 

this  is  the  same  approximation  that  underlies  Eq.  (3.21) .  By  using  z  x  p  =  $  and  cos  p0  p+ sin  p0z  = 
ft,  we  obtain  then  the  conventional  UTD  result 

Hd~Hd(UTD)  =  -Uifca  {DsHo-W  +  D^PPo  H1}  ,  (3.29) 

analogous  to  Eq.(3.15). 

To  obtain  the  exact  expression  for  the  diffracted  field,  we  differentiate  the  diffraction  coefficients 
in  Eq.(3.25)  using  the  relation  (3.8),  and  express  the  result  in  terms  of  the  unit  vectors  sd  =  s,  d>, 
and  p  by  writing  p  =  sin  p0s  +  cos  p0  p.  We  find  that  exact  diffracted  field  is 

Hd  =  Hd(UTD)  +  <5Hd  f  (3.30) 

where  the  additional  term  <5Hd  can  be  represented  in  terms  of  a  new  dyadic  diffraction  coefficient 

V  as  -i 

SHd  =  4=eifcsX»Hi  .  (3-31) 

Vs 

Here  V  is  a  dyadic  form  which  can  be  expressed  in  terms  of  the  unit  vectors  associated  diffracted 
and  incident  rays.  In  a  compact  matrix  notation  it  has  the  form 

^ P<p  Vpp 

'D<t><p  _  A 

(3.32) 

where  the  dyadic  coefficients  associated  with  the  vectors  p  and  4>  are 


Vp4>  =  - C(k,P0 )  {sin V*—  [l  -  F  (£_)]  -  sin^+  [l  -F  ($+)] }  ,  (3-33a) 

=  -C(k,p0)  sin2p0  {cos^_  [l  -  F  (£_)]  -  cos^+  [l-F  (£+)] }  ,  (3.33b) 

Vpp  =  -C(k,P0)  sin  2P0  cos  P0  {cosip_  [l-F  (£_)]  +  cos  ii>+  [l-F  (£+)] }  ,  (3.33c) 

=  C(k,P0)  cosP0  {sin  if>_  [l-F  (£_)]  +  sin  if>+  [l-F  (£+)] }  .  (3.33d) 


A  r,  n  sinA>  0 

<p  o 

v  =  [s  4>  0)  =  [*  4>  0]  0  1 

Vpt  Vpp.  L  °J  [cos A)  0 


Here 


ei7r/4 

C(k,p0)=  — - 

V  27rfc  sin  p0 


(3.34) 


^  =  I(^  T  ^0),  and  ^  =  2k s  sin 2p0  cos2  We  note  that  the  diffracted  field  <5Hd  is  not 

transverse  -  it  contains  a  longitudinal  component  proportional  to  s.  We  also  note  that  all  the 
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diffraction  coefficients  (3.33)  are  proportional  to  [1  -  F(-  ■  •)]>  i.e.,  they  fall  off  rapidly  outside  the 
transition  region. 

In  the  special  case  of  normal  incidence  (/ 30  =  7r/2)  the  total  diffracted  field  can  be  written 
explicitly  as 

Hd(p,  4>)  =  -e1,r/4eifcp  yL=r=  {P  [2  sinV’.  [l  -  F  (£_)]  -  2  sin^  [l  -F  (£+)]] 

*  P+  4>  [2  cosip_  [1  -  F  (£_)]  -  2  cos  ^4.  [1  -  F  (£+)] 

+  sec  ip_  F(£_)  -  sec  ip  +  F(£+)] 

+  z  [s ecip_  F(£_) +  secip+F(£+)]  Hp]  ,  (3.35) 


where  now  £,  =  2  kp  cos2  ip=.  The  magnetic  field  discontinuity  is  then 


AHd(p)  =  -ei,r/4eifcp 


V27T  Ttp 


j-p  2  sin  —  1  -  F  ^2  kp  cos2  +  z  sec  yF  ^2  kp  cos2  J  • 


(3.36) 


Here  the  second  term,  proportional  to  z,  is  the  same  as  in  conventional  UTD,  and  is  due  to  “hard 
diffraction”.  It  represents  the  magnetic  field  directed  along  the  edge,  giving  rise  to  the  current 
(3.10)  in  the  ^-direction, 

JM  =  4utdKp)  =  -2AHM,  <3-37> 

normal  to  the  edge.  The  first  term,  proportional  to  p  =  x  is  the  new  near  field  contribution  due  to 
our  modification  of  UTD.  It  is  associated  with  “soft  diffraction”,  and  yields  a  current  component 

Jz(p)  =  6Jz(p)  =  2AHp(p)  (3-38) 

in  the  ^-direction,  i.e.,  along  the  edge.  As  expected,  this  current  component  behaves  for  small  p  as 
1/y/p,  and  decreases  rapidly  away  from  the  edge  (outside  the  transition  region). 

In  Fig.  19  we  plot  real  parts  of  the  current  components  Jx  and  Jz  as  functions  of  kp,  for  two 
values  of  <p0:  90°  and  170°.  The  plots  show  that,  as  <p0  approaches  180°  the  transition  region 
extends  to  larger  values  of  kp,  and  the  current  Jz  remains  sizable  over  a  wider  range. 

It  is  of  interest  to  note  that  the  size  of  our  correction  <5Hd  to  the  diffracted  field  is  of  the  same 
order  as  the  difference  between  the  UTD  and  GTD  result  (the  latter  being  UTD  with  all  transition 
functions  F  replaced  by  unity),  i.e.,  schematically, 

<5Hd  ~  Hd  (UTD)  -  Hd  (GTD)  .  (3-39) 

Indeed,  both  quantities  are  sums  of  terms  proportional  to  [1  —  F{-  ■  •)].  For  instance,  in  our  example 
Eqs.  (3.37)  and  (3.38)  imply  that,  assuming  unit  values  of  Hp  and  H^, 

SJM  =  2  sin^o  [4'™)W  -  4GTD,W]  ■  <3-4°> 
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For  a  more  quantitative  estimate  we  evaluate  the  integrals  of  the  current  components,  e.g., 


The  function  g(a)  is  defined  for  -oo  <  a  <  \  with  the  square-root  branch  specified  by  Ima  <  0, 
i.e.,  y/l  -  =  -i  -  1  f°r  a  <  -1.  It  is  real,  and  is  normalized  so  that  5(0)  =  \  (corresponding 

to  <f>0  —  7r/2),  and  g(-oo)  =  1  (corresponding  to  <f>0  =  tt).  From  its  plot  (Fig.  20)  it  is  apparent 
that  the  size  of  both  SJZ  and  4UTD)  -  «/iGTD)  increases  when  we  approach  <j>0  =  180°. 


Fig.  19:  Real  parts  of  the  current  com-  Fig.  20:  The  factor  g  of  Eq.(3.42)  plot- 

ponents  (3.37)  and  (3.38)  as  functions  of  ted  as  a  function  of  the  angle  <j>Q. 

kp,  for  <f>0  =  90°  and  <f>0  =  170°,  assum¬ 
ing  unit  values  of  ifg  and  H^. 

In  Section  4  we  present  examples  in  which  we  compare  scattering  cross-sections  obtained  using 
the  MoM,  with  those  computed  using  currents  determined  from  GTD,  UTD,  and  our  modified 
version  of  UTD.  We  find  there  the  contributions  due  to  the  modification  of  UTD  significantly 
improve  the  agreement  of  asymptotic  approximations  with  the  exact  results. 
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3.2.  Direct  field  computation 

We  return  now  to  the  problem  of  determining  scattered  fields  directly  form  the  fields  defined  on 
the  WFs. 

Properties  of  the  WFs  reflected  from  smooth  curved  surfaces,  and  surfaces  consisting  of  flat 
facets,  are  quite  different.  In  the  first  case  the  WF  is  curved,  and  the  individual  rays  have  also 
nonzero  curvatures.  In  the  second  case,  meanwhile,  the  reflected  WF  consists  of  flat  surface 
segments,  and  the  rays  carry  zero  curvatures  (provided  we  consider  a  plane  incident  WF). 

Similarly,  WFS  diffracted  from  straight  and  curved  edges  behave  differently:  they  are  either 
cylindrical  (conical),  or  spherical  waves.  In  the  following  we  consider  thus  only  the  processes  which 
generate  spherical  waves,  amenable  to  direct  computation  of  the  scattering  cross-section. 

By  definition,  the  scattering  amplitude  is  defined  as  the  appropriately  normalized  asymptotic 
value  of  the  scattered  field.  A  vector-valued  scattering  amplitude  can  be  defined  in  terms  of  the 
magnetic  fields  as  the  limit 

A(R)=  lim  47ri?e-ifc*H(R)  ,  (3.43) 

v  R-+  oo 

assuming  a  unit  strength  incident  field,  |Hin|  =  1,  and  using  the  field  H  including  the  rapidly 
varying  phase  factor  (as  opposed  to  the  smoothly  varying  field  H  used  in  Eq.(3.2)).  According 
to  the  laws  of  GO  we  can  calculate  the  field  H  at  any  observation  point  R  provided  we  identify 
all  rays  passing  through  that  point.  Fig.  21  shows  two  such  rays,  originating  at  points  rx  and 
r2  of  a  certain  WF,  characterized  by  the  given  fixed  phase  ks  (where  s  denotes  the  length  of  the 
propagation  path). 


Fig.  21:  Contributions  to  the  field  at  the  observation  point  R  due  to  two  rays  originating  at 
points  Tj  on  segments  Xj  of  a  WF  {j  =  1,2),  and  propagating  along  the  normals  n*  towards  the 
observation  point. 


Since  the  rays  in  the  regions  outside  the  scatterer  evolve  according  to  the  laws  of  GO,  the 
contributions  of  the  fields  at  points  Tj  to  the  field  H(R)  are  given  by  (cf.  Eq.(3.2)) 


h(j)(R) 


JJ)  Ji) 
Pi  P2 


\  (Pl]  +  Sj)  (P2  +  Sj) 


eiks*  H(r,)  , 


(3.44) 
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where  p^  are  the  ray’s  curvatures  at  the  points  rj5  and  Sj  =  |R  -  r^l  is  the  propagation  path 
length.  Since,  in  the  limit  R  ->  oo,  we  have  Sj-R-  Rrp  the  amplitude  (3.43)  becomes 

«4(R)  =  4t r  y/p?  p(2j)e'ikA^  H^)  ;  (3-45) 

j 

actually,  in  practice  we  apply  the  GO  evolution  laws  to  real  fields  with  the  common  phase  exp(ifcs) 
factored  out. 

The  conventionally  normalized  differential  scattering  cross-section  (per  unit  solid  angle)  for 
scattering  in  the  direction  R  and  for  the  scattered  wave  magnetic-field  polarization  h  is  then  given 

by  1 

<re(R)  =  ^|hVl(R)|2.  (3-46) 

According  to  Eq.(3.45),  the  scattering  amplitude  has  the  correct  dimension  of  length,  and  the 
scattering  cross-section  the  dimension  of  length  squared. 

The  procedure  for  computing  the  scattering  amplitude  -4(R)  is  thus  as  follows: 

1.  Select  a  sufficiently  evolved  WF  characterized  by  some  evolution  parameter  s.  ^ 

2.  On  this  WF  locate  all  points  rj  whose  normals  coincide  with  the  direction  R. 

3.  Sum  contributions  of  these  points  according  to  Eq.(3.45). 

Having  found  the  vector-valued  amplitude  A( R),  we  can  then  computing  the  conventional  scat¬ 
tering  amplitudes  by  taking  scalar  products  of  A  with  the  magnetic  field  polarization  vectors  h  of 
the  scattered  wave.  Consequently,  the  differential  polarized  cross-sections  are  proportional  to  the 
squares  of  the  absolute  values  of  such  amplitudes. 

Some  remarks  are  in  order  here: 

(A)  In  the  step  1.  above  by  a  “sufficiently  evolved”  WF  we  mean  a  WF  containing  contributions 
from  a  sufficient  number  of  reflections  (if  there  are  multiple  reflections  in  the  problem).  In 
practice,  we  can  compute  the  scattering  amplitudes  for  a  number  of  WFs,  and  monitor  the 
convergence  of  the  result.  Typically,  the  convergence  is  relatively  fast,  since  the  consecutive 
reflections  give  weaker  and  weaker  contributions  to  the  field. 

(B)  Since  the  WF  is  represented  in  terms  of  a  finite  number  of  rays,  we  cannot  expect  to  find 
normals  n^  pointing  exactly  in  the  direction  R.  Instead,  we  have  to  use  an  interpolation 
procedure:  We  first  map  the  mesh  defining  the  WF  onto  a  unit  sphere,  on  which  the  nodes  are 
the  normals  n.  On  this  unit  sphere  we  find  then  all  the  mesh  faces  (triangles)  intersected  by 
the  direction  vector  R.  Finally,  for  each  such  triangle  we  interpolate  the  field  values  associated 
with  the  three  vertices  to  the  point  representing  the  required  direction  R;  for  sufficiently  small 
triangles,  linear  interpolation  is  adequate.  An  important  element  of  the  procedure  is  that  in 
the  interpolation  we  use  the  smoothly  varying  fields  H  (as  in  Eq.(3.2)),  and  only  then  we 
multiply  the  result  by  the  appropriate  phase  factor. 
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3.3  Surface  current  computation 

The  goal  of  the  current-based  approach  is  to  obtain  equivalent  surface  currents  corresponding  to 
the  field  of  the  WF,  and  use  these  currents  to  compute  the  scattered  field,  or  for  other  purposes 
(as  we  discuss  below).  For  a  p.e.c.  scatterer  the  equivalent  surface  electric  current  is  given  by 
Eq.(3.1),  hence  the  values  of  the  magnetic  field  on  the  surface  are  required. 

Besides  computing  the  scattered  field  (as  the  Kirchhoff  integral  of  the  current),  a  natural 
applications  of  this  type  of  calculations  might  be  to  develop  an  interface  to  integral-equation 
methods,  in  which  the  unknown  quantities  are  surface  currents.  There  is  a  need  for  such  an 
interface  in  several  circumstances: 

1.  Hybrid  methods  require  input  for  solving  low-frequency  integral  equations  in  the  form  of 
the  initial  current  generated  by  the  field  originating  from  the  high-frequency  components  of 
the  object.  Computation  of  this  current  provides,  effectively,  a  coupling  between  the  system 
components  described  by  low-frequency  and  asymptotic  methods. 

2.  In  an  “enhanced”  formulation  of  Physical  Optics  (PO)  the  required  input  for  scatterer  field 
computation  is  the  current  induced  on  the  scatterer  surface  by  the  incident  field,  and  by  the 
fields  due  to  GO  reflections  (in  this  case  the  WF  is  thus  obtained  in  the  GO  approximation). 

3.  The  rigorous  method  of  high-frequency  asymptotic  integral  equations  considered  in  [2, 3]  re¬ 
quires  an  Ansatz  for  the  solution  for  the  surface  current.  This  Ansatz  is  in  the  form  of  a 
linear  combination  of  unknown  slowly  varying  functions  multiplying  rapidly  oscillating  func¬ 
tions  representing  asymptotic  high-frequency  solutions.  These  asymptotic  solutions,  generated 
by  all  possible  reflection  and  diffraction  processes,  can  be  very  naturally  extracted  from  the 
purely  geometrical  construction  of  WFs  (even  without  the  need  of  computing  diffracted  field 
by  means  of  UTD  or  similar  methods). 

One  of  the  first  questions  we  encounter  when  trying  to  determine  surface  currents  is  how  these 
currents  should  be  represented  and  parameterized.  At  high  frequencies  it  is,  clearly,  impractical  to 
use  a  representation  of  the  current  in  terms  of  their  values  at  selected  points  on  the  surface,  since 
this  would  require  about  10  sampling  points  per  wavelength.  A  more  useful  representation  is  its 
parameterization  as  a  sum  of  rapidly  oscillating  exponential  factors  multiplied  by  smoothly  varying 
coefficients.  Such  parameterizations  would  be  applicable  in  some  surface  patches  Dj  centered 
around  some  “current  points”  Rj.  The  sizes  of  the  regions  Dj  is  assumed  to  be  of  the  order  of  the 
local  ray  spacing,  and  they  are  required  to  cover  the  entire  scatterer  surface  S. 

Let  us  assume  for  a  moment  that  a  set  of  M  distinct  rays  (originating  from  different  segments 
of  the  WF)  coalesce  at  a  current  point  Rj.  The  conventional  WF  form  of  the  fields  suggest  then 
the  representation 

M 

J(r)  =  Am(r)eifcSm(r)  for  r  G  Dj  >  (3-47) 

m=l 

i.e.,  a  sum  of  contributions  of  distinct  rays  arriving  at  R,.  Here  5m’s  are  the  phases  of  the 
individual  rays,  defined  so  that  VSm(r)  =  nm,  where  nm  is  the  unit  vector  in  the  direction  of  the 
ray  propagation.  The  amplitudes  Am  are  parameterized  as  constant  amplitudes  multiplied  by  the 
conventional  ray  divergence  factors  involving  the  ray  evolution  parameter  and  the  curvatures. 
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If  the  rays  arrive  at  the  surface  not  exactly  at  the  “current  point”  Rj,  but  in  its  vicinity, 
it  is  necessary  to  devise  a  procedure  to  interpolate  the  rays  to  the  point  Rj,  and  then  find  a 
parameterization  of  the  phases  and  coefficients  in  Eq.(3.47).  Such  a  procedure  can  be  based  on 
WF  field  interpolation  within  sets  of  neighboring  rays,  similarly  to  methods  using  “ray  cells”  (i.e., 
prisms  formed  by  the  corresponding  faces  of  the  meshes  on  the  two  consecutive  WFs). 

Our  algorithm  can  be  formulated  as  a  procedure  for  computing  field  values  at  any  set  of 
observation  points  (not  necessarily  points  on  the  scatterer  surface).  Currents  can  be  then  obtained 
by  applying  the  algorithm  to  surface  points,  and  using  J  =  2n  x  AH,  where  AH  is  the  discontinuity 

of  the  magnetic  field  across  the  surface. 

Conceptually,  the  interpolation  algorithm  can  be  stated  as  follows: 

1.  Specify  a  set  K  of  “field  points”  Rj  at  which  the  fields  are  to  be  computed;  the  spacing  of 
these  points  should  be  sufficiently  small  to  resolve  the  variation  of  the  amplitudes  of  the  fields 
(but  not  variation  of  their  phases). 

2.  At  every  step  of  the  WF  evolution  find,  for  each  “field  point”  Rj,  a  set  of  nearby  rays  (within 
the  distance  of  a  few  average  ray-ray  spacings).  These  rays  are  selected  from  the  set  of  rays 
propagating  from  the  previous  WF  to  the  current  WF  x  (Fig-  22)- 

3.  For  each  of  nearby  rays  £  construct  a  ray  £R  passing  through  the  point  Rj  and  parallel  to 

j 

the  ray  £.  The  length  As  of  the  ray  is  set  to  the  distance  between  the  point  R,  and  the  WF, 
measured  along  the  direction  of  the  ray  £. 

4.  By  evolving  the  ray  £R  forward  in  time  find  its  intersection  with  the  WF  %.  Identify  the 
WF  face  /  intersected  by  the  ray  (Fig.  22).  (We  note  that  a  continuation  of  the  ray  £R_  may 
intersect  further  segments  of  the  WF,  but  theses  intersections  are  irrelevant,  since  we  are  only 
computing  contributions  of  the  rays  evolved  up  to  the  WF  %•) 

5.  Identify  the  rays  £fv  £/2,  and  £/3,  associated  with  the  corners  of  the  face  /.  Evolve  these  rays 
backward  through  the  distance  As  to  form  an  image  /'  of  the  face  /.  We  refer  to  the  prism 
built  on  the  faces  /  and  /'  as  a  ray  tube. 

6.  Check  if  the  point  Rj  is  located  inside  the  face  If  it  is,  continue  to  the  point  7.  If  it  is 
not,  determine  the  average  direction  n123  of  the  rays  £fi  (i  =  1, 2, 3),  construct  a  new  ray  £R^ 
emerging  from  Rj  in  the  direction  of  ^ 1 2 3  ■  &fld  return  to  the  point  4. 

7.  Interpolate  field  values  associated  with  the  ends  of  the  rays  (fi  to  the  point  Rj. 
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Fig.  22:  The  algorithm  for  computing  a  contribution  of  a  ray  to  the  fields  at  the  observation  point 

Rj,  by  interpolating  between  rays  an<^  £/3  associated  with  a  face  /  on  the  WF  x •  The 

heavy  arrow  denotes  the  direction  of  WF  propagation. 

Several  remarks  and  clarifications  are  in  order  here: 

(a)  The  steps  3  to  7  are  repeated  for  all  rays  passing  near  the  observation  point  Rj.  If  several 
near  rays  yield  the  same  ray  tube,  only  one  of  them  is  taken  into  account.  Contributions  to 
the  fields  from  all  distinct  ray  tubes  are  added  to  the  field  representation  at  Rj.  The  field  is 
represented  in  analogy  to  Eq.(3.47),  i.e.,  in  the  form  of  constant  vector  amplitudes  and  the 
corresponding  wave  vectors.  Specifically,  for  each  of  the  ray  tubes  containing  the  point  R 
we  store  the  field  amplitudes  Em(Rj)  and  the  phase  Sm( Rj)  evaluated  at  the  observation 
point,  as  well  as  the  gradient  njiTn  =  VSm(Rj)  at  that  point.  These  data  provide  a  local 
parameterization  of  the  field  in  the  neighborhood  of  the  point  Rj  as  a  sum  of  plane  waves. 

(b)  If  the  observation  points  R^  are  located  on  the  scattered  surface,  it  is  understood  that  the 
WF  x  is  the  WF  computed  as  in  the  absence  of  the  scatterer,  i.e.,  before  the  effect  of  rays’ 
reflections  is  taken  into  account.  In  fact,  in  our  present  implementation  the  WFs  are  first 
evolved  without  reflections,  and  then  transformed  by  reflecting  the  rays. 

(c)  We  note  that  our  algorithm  for  field  interpolation  differs  from  the  method  being  used  in 
previous  WF  approaches,  and  based  on  a  somewhat  different  concept  of  a  i4ray  cell  [5, 6].  In 
that  construction  of  a  ray  cell  is  specified  in  terms  of  two  consecutive  WFs  (the  “previous” 
and  the  “current”  ones),  while  our  algorithm  employs  only  one  WF  (the  “current”  WF).  Our 
procedure  was  mainly  motivated  by  applications  to  diffraction:  in  this  case  we  may  encounter 
a  situation  where  the  “previous”  WF  simply  does  not  exist,  because  the  current  WF,  due  to 
a  diffraction  process,  was  newly  created  in  the  considered  evolution  step. 

(d)  In  order  to  compute  the  induced  electric  field  on  an  open  surface,  it  is  necessary  to  evaluate 
the  field  slightly  above  and  slightly  below  the  surface.  In  doing  that,  we  have  to  check  which 
side  of  the  surface  is  illuminated  by  the  individual  rays  of  the  WF,  and  we  have  to  take  into 
account  possible  discontinuities  of  diffracted  waves  across  the  surface. 
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(e)  In  the  interpolation  of  fields  associated  with  the  rays  we  have  to  take  into  account  the  fact 
that  the  fields  may  not  be  ray-optical  (in  the  UTD  transition  regions,  including  the  region 
near  the  diffraction  source). 

To  summarize,  the  algorithm  for  computing  the  surface  current,  as  described  above,  is  similar  to 
the  previous  method  for  computing  the  asymptotic  far  field,  except  that  now  the  observation  points 
are  located  on  the  scatterer  surface.  The  resulting  complication  is  that  it  is  necessary  to  take  into 
account  contributions  of  all  WFs  during  their  evolution,  rather  than  only  the  far  evolved  WFs, 
after  all  relevant  scattering  processes  have  occurred.  In  addition,  an  efficient  current  representation 
requires  local  analytic  parameterizations  of  the  current,  rather  than  merely  field  or  current  values 
at  some  sampling  points. 

We  also  note  that  the  above  current  computation  algorithm  is  applicable  both  to  WFs  resulting 
from  reflection  as  well  as  diffraction  effects.  It  can  be  thus  used  both  to  generate  currents  as  input 
for  Kirchhoff  integral-type  computations,  and  for  creating  the  initial  Ansatz  (including  diffraction) 
for  the  solution  of  rigorous  asymptotic  surface-current  integral  equations. 

3.4.  Assessment  of  the  scattered  field  computation  procedures 

From  the  practical  point  of  view,  the  main  advantage  of  the  direct  computation  of  the  scattered 
fields  is  the  relative  simplicity  of  the  algorithm.  Its  drawbacks,  however,  are  the  difficulties  which 
appear  near  the  scattered  field  caustics  (we  will  see  an  example  of  that  for  scattering  on  a  circular 
disc),  and  the  fact  that  the  only  contribution  to  the  cross-section  comes  from  processes  generating 
asymptotically  spherical  waves .  Thus,  for  example,  in  scattering  on  a  polyhedral  object  with  flat 
faces,  corner  diffraction  is  the  only  mechanism  generating  asymptotically  spherical  waves,  hence, 
corner  diffraction  is  the  only  process  contributing  to  the  scattered  field  and  to  the  cross  section 
(edge  diffraction  generates  cylindrical  or  conical  waves).  This  situation  appears  rather  unnatural, 
considering  the  fact  that  edge  diffraction  is  certainly  an  important  scattering  process. 

The  advantage  of  the  procedure  using  induced  currents  to  compute  the  scattered  fields  is  its 
applicability  to  all  diffraction  processes,  and  the  fact  that  it  avoids  difficulties  due  to  caustics  in  the 
scattered  field.  (The  only  difficulty  remains  if  the  caustic  is  located  on  or  intersects  the  scatterer’s 
surface.)  A  certain  disadvantage  of  the  approach  is  a  relative  complexity  of  the  algorithm  for 
current  computation. 

Below  we  present  results  obtained  using  both  procedures  for  computing  the  scattered  fields. 


4.  Examples 

We  give  now  results  of  scattered  field  and  cross-section  computations  for  several  simple  scatterers 
of  various  types. 

Reflection  off  a  system  of  spheres 

We  consider  multiple  reflection  process  on  a  system  consisting  of  two  spheres  of  unit  radii,  separated 
by  a  distance  equal  to  the  sphere  diameter,  as  shown  schematically  in  Fig.  23.  The  incident  plane 
WF  approaches  the  spheres  along  the  negative  2- axis,  with  the  electric  field  along  the  2;- axis. 
The  evolution  step  length  is  0.3.  For  the  cross-section  computation  we  assume  the  wavelength 
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A  =  0.2  (i.e.,  the  spheres  diameters  10  A).  Fig.  24  shows  the  differential  bistatic  cross-section  for 
the  vertical  ( 6 -)  polarization,  compared  with  the  exact  result  obtained  using  a  LF  code  (in  this 
case,  our  integral-equation  solver  with  FMM  compression),  with  the  discretization  leading  to  a 
problem  with  about  100,000  unknowns.  The  differential  cross-section  shown  is  the  conventional 
dimensionless  cross-section  &  ==  of  A2  normalized  to  the  wavelength  squared. 

The  rapidly  oscillating  structure  in  the  cross-sections  of  Fig.  24  is  due  to  interference  of 
reflections  from  the  two  spheres.  The  flat  segment  in  the  cross-section  results  from  one  sphere 
obscuring  the  view  of  the  other,  and  is  due  to  single  reflection  from  one  sphere  only.  The  agreement 
of  the  asymptotic  (GO)  and  exact  cross-sections  is  good  in  the  near-back-scattering  region,  where 
the  main  contributions  come  from  single  reflections  from  the  individual  scatterers.  As  expected, 
the  cross-sections  do  not  agree  near  90°  and  180°,  where  diffraction  effects  are  important. 


Fig.  23:  A  multiple-reflection  problem 
with  two  spheres. 


scattering  angle 


Fig.  24:  Comparison  of  the  asymptotic 
(GO)  cross-section  with  the  exact  LF 
code  result. 


Diffraction  on  a  circular  disc 

We  consider  here  scattering  on  a  circular  disc  of  radius  o,  placed  in  the  (x,y)  plane,  illuminated 
from  above  with  an  incident  plane  wave,  with  the  electric  field  along  the  x  axis.  For  simplicity,  we 
take  the  normal  incidence;  in  this  case  the  first-order  diffraction  contribution  to  the  GTD  scattering 
amplitude  in  the  scattering  plane  <f>  =  0  is  [13, 14] 

cos(fcq  sin 6  -  7t/4)  .  sin(fca  sing  -  7r/4)  /4 

cos  6/2  sin  6/2 

This  expression  is  valid  away  from  the  specular  reflection  direction  6  =  0  (which  is  also  a  caustic  of 
the  reflected  and  diffracted  rays).  Eq.(4.1)  is  due  to  the  interference  of  two  rays  diffracted  from  the 
opposite  points  on  the  disc  edge,  both  located  on  the  x  axis.  We  also  note  that  for  computation  of 
the  far-field,  away  from  the  caustics  and  transition  regions,  GTD  and  UTD  give  identical  results. 


A(6)  =  2]j 


2  na 
k  sin# 
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In  Fig.  25  we  compare  the  GTD/UTD  cross-section  based  on  Eq.(4.1)  with  the  result  of  our 
WF  computation  (according  to  Eqs.  (3.45)  and  (3.46)),  and  with  the  rigorous  MoM  result.  The 
computations  were  done  for  a  disc  of  radius  a  =  5  A. 


Fig.  25:  The  differential  bistatic  cross-section  (in  dB,  normalized  to  the  wavelength  squared)  for 
a  circular  disc,  computed  using  GTD  (Eq.(4.1)),  the  Method  of  Moments  (MoM),  Physical  Optics 
(PO),  and  our  WF  evolution  algorithm. 

The  comparison  shows  that  the  agreement  of  the  WF  calculation  with  GTD  (away  from 
the  backscattering  direction)  is  very  good,  except  near  the  first  (very  narrow)  minimum.  The 
discrepancy  in  this  case  is  likely  due  to  inaccuracies  in  field  interpolation  in  the  range  where  the 
field  varies  very  rapidly. 

As  expected,  neither  GTD/UTD  nor  WF  computation  reproduce  correctly  the  main  diffrac¬ 
tion  peak  at  0  =  0  (which  occurs  in  the  direction  of  the  reflected  and  diffracted  rays  caustic). 
Clearly,  GTD/UTD  give  an  infinite  cross-section  at  0  =  0.  This  deficiency  may  be  corrected  by 
our  alternative  approach  to  scattered  field  computation  based  on  evaluating  the  surface  current 
distribution  (provided  the  current  itself  can  be  reliably  calculated  in  the  caustic  region. 

At  the  same  time,  the  conventional  PO  approximation  is  expected  to  fail  at  larger  scattering 
angles  (in  the  bistatic  cross-section).  This  well  known  fact  is  confirmed  by  the  plots  shown  in 
Fig.  25.  Thus,  away  from  the  specular  reflection  and  caustic  directions,  GTD/UTD  gives  predic¬ 
tions  definitely  better  than  PO  —  unless  the  latter  approximation  is  improved  by  adding  diffr action 
contributions  to  the  current  itself. 

The  discrepancies  between  the  exact  (MoM)  result  and  GTD/UTD  or  WF  calculations  are 
evidently  due  to  the  fact  that  latter  include  only  the  first-order  diffraction .  Repeated  multiple 
diffraction  from  the  edges  of  the  same  flat  screen  can  be  included  in  the  WF  computation,  but 
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requires  &n  extension  of  the  diffraction  &lgonthm  described  above.  The  present  nlgorithm  is  not 
able  to  detect  an  edge  as  giving  rise  to  diffraction  if  that  edge  is  illuminated  by  a  ray  strictly 
tangential  to  the  screen;  this  problem  is  one  of  the  issues  we  discuss  in  Section  5. 

Diffraction  on  a  rectangular  plate 

In  this  example  we  consider  a  6 A  x  6 A  square  plate  placed  in  the  {x,y)  plane,  illuminated  by  a 
plane  wave  of  vertical  (0)  or  horizontal  (</>)  polarization,  with  the  incident  wave  vector  in  the  (a;,  z) 
plane  (<f>  =  0).  Here  we  evaluate  the  scattered  field  by  first  computing  the  induced  current  on  the 
scatterer  surface,  using  the  algorithm  described  in  Section  3.3.  In  evaluating  the  current  we  use 
either  only  the  incident  wave  (is  equivalent  to  the  PO  approximation),  or  include  contributions 
from  single  edge  diffraction  utilizing  GTD,  UTD,  and  “near-field  UTD”  (nUTD),  as  described  in 
Section  3. 

Figs.  26  and  27  show  the  results  for  the  vertical-polarization  bistatic  cross-section  as  the 
function  of  the  scattering  angle  0,  for  two  incidence  angles:  normal  incidence  (05  =  0°)  and  near¬ 
grazing  incidence  (0i  =  80°). 

In  the  normal  incidence  case  for  the  vertical  polarization  (Fig.  26a)  the  cross-section  resulting 
from  the  PO  current  is  much  too  low  for  theta  approaching  90°.  The  GTD  current  yields  a  better 
result,  while  the  UTD  and  nUTD  currents  give  rise  to  cross-sections  almost  indistinguishable  from 
the  MoM.  For  the  horizontal  polarization  (Fig.  26b),  on  the  other  hand,  PO,  GTD,  and  UTD 
currents  yield  almost  identical  cross-sections,  which,  however,  do  not  agree  with  the  MoM  result 
neat  0  =  90°.  The  modified  UTD,  nUTD,  provides  a  much  better  agreement  with  the  MoM. 
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Fig.  26a:  The  differential  bistatic  cross 
section  (in  dB,  normalized  to  the  wave¬ 
length  squared)  for  the  6  A  x  6  A  plate 
at  normal  incidence  (0*  =  0°)  and  for 
vertical  polarization,  computed  by  us¬ 
ing  the  MoM  and  our  current-based  al¬ 
gorithm  utilizing  PO,  GTD,  UTD,  and 
nUTD  currents. 


Fig.  26b:  Same  as  Fig.  26a,  but  for  hor¬ 
izontal  polarization. 


In  the  near-grazing  incidence  case  for  the  vertical  polarization  (Fig.  27a)  the  PO  current  gives 
a  poor  agreement  with  the  MoM  for  all  scattering  angles.  The  GTD  current  results  in  an  even 
worse  approximation.  The  remaining  procedures  (UTD  and  nUTD)  provide,  on  the  other  hand, 
very  good  agreement  with  the  MoM.  In  the  horizontal  polarization  case  (Fig.  27b)  the  PO,  GTD, 
and  UTD  currents  give  fairly  close  results,  but  all  strongly  disagree  with  the  MoM.  The  nUTD 
provides  a  much  better  agreement. 

We  also  note  a  very  poor  agreement  of  the  GTD  result  for  vertical  polarization  (Fig.  27a), 
due  to  the  incorrect  behavior  of  the  GTD  fields  in  the  transition  region  (which,  for  near-grazing 
incidence,  extends  over  a  large  area  of  the  plate) . 
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Fig.  27a:  The  differential  bistatic  cross 
section  for  the  6  A  x  6  A  plate  at  near¬ 
grazing  incidence  (0i  —  80°)  and  for 
vertical  polarization,  computed  by  us¬ 
ing  the  MoM  and  our  current-based  al¬ 
gorithm  utilizing  PO,  GTD,  UTD,  and 
nUTD  currents. 


Fig.  27b:  Same  as  Fig.  27a,  but  for  hor¬ 
izontal  polarization. 


Results  for  the  back-scattering  cross-section  are  shown  in  Fig.  28.  The  features  of  the  various 
approximations  are  here  similar  to  those  seen  for  bistatic  cross-sections.  We  note,  in  addition,  the 
discontinuity  of  the  GTD  vertical-polarization  cross-section  near  6  =  90°,  due  to  the  nonuniform 
behavior  of  the  GTD  diffraction  coefficients  in  the  transition  region. 
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Fig.  28a:  The  differential  back-scat-  Fig.  28b:  Same  as  Fig.  28a,  but  for  hor- 

tering  cross  section  for  the  6  A  x  6  X  plate,  izontal  polarization. 

for  vertical  polarization,  computed  by 

using  the  MoM  and  our  current-based 

algorithm  utilizing  PO,  GTD,  UTD,  and 

nUTD  currents. 

As  a  further  illustration  we  show  in  Fig.  29a  the  distribution  of  the  induced  current  Re  Jx 
(computed  using  nUTD)  for  the  case  of  the  near-grazing  incidence  of  Fig.  27b,  with  the  v-polarized 
incident  wave. 


Fig.  29a:  Distribution  of  the  current  Fig.  29b:  The  same  as  Fig.  29a,  but  for 

component  Re  Jx  for  a  v-polarized  wave  an  h-polarized  incident  wave, 

at  the  incidence  angle  6i  =  80°. 

The  current  Jx  consists  of  the  incident  wave  (PO)  and  diffraction  contributions;  the  latter  is 
mainly  due  to  diffraction  on  the  “trailing  edge” ,  and  results  in  the  enhancement  (compared  to  PO) 
of  the  back-scattering  cross-section  at  near-grazing  angles,  as  seen  in  Fig.  28.  (That  enhancement 
is  commonly  termed  the  “surface  wave”  contribution.) 
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For  comparison,  Fig.  29b  shows  the  distribution  of  the  same  current  component  as  Fig.  29a, 
but  for  an  h-  (<£-)  polarized  incident  wave.  In  this  case  the  contributions  to  R eJx  come  exclusively 
from  diffraction  on  the  edges  parallel  to  the  x-axis;  it  can  be  seen  that  the  currents  are  nonzero 
only  inside  the  relevant  diffraction  cones. 

The  features  of  the  approximate  cross-sections  in  Figs.  26  and  27  can  also  be  understood  by 
analyzing  the  behavior  of  the  induced  currents.  In  particular,  the  large  difference  between  the 
nUTD  and  the  remaining  approximate  results  in  Figs.  26b  and  27b  is  due  to  the  fact  that  only 
nUTD  correctly  reproduces  the  behavior  of  the  currents  flowing  in  the  y-direction  near  the  edges 
aligned  with  the  y-axis. 

We  have  analyzed  many  other  cases  of  scattering  on  plates  (for  other  cross-section  cuts,  and 
other  incidence  angles)  and  we  found  that  the  currents  computed  using  nUTD  consistently  provide 
an  improvement  relative  to  both  GTD  and  conventional  UTD.  Nevertheless,  the  agreement  with 
the  rigorous  results  is  often  not  as  good  sis  that  shown  in  Figs.  26-28. 

There  are  many  possible  sources  of  disagreement  between  the  MoM  currents  and  currents 
evaluated  using  an  asymptotic  theory.  We  stress  here  that  nUTD  applied  to  a  rectangular  plate 
reproduces  exactly  the  solutions  of  canonical  (half-plane)  problems  for  the  four  plate  edges.  There¬ 
fore,  the  most  likely  missing  diffraction  mechanism  is  corner  diffraction,  understood  not  only  as 
a  term  in  the  asymptotic  expansion  of  the  Kirchhoff  integral,  but  also  as  due  to  specific  physical 
diffraction  processes  occurring  in  the  presence  of  a  corner.  Unfortunately,  the  known  solutions  of 
the  canonical  problem  (a  plane  angular  sector)  are  only  given  as  a  rather  complicated  series  of  Lame 
functions  in  the  sphero-conal  coordinate  system  [17, 18, 19].  For  this  reason,  the  comer  diffraction 
problem  has  not  yet  been  completely  solved,  although  asymptotic  forms  of  the  corner-diffracted 
field  (corner  diffraction  coefficients)  have  been  obtained  under  some  conjectures  about  the  solution 
[20, 21],  or  for  the  rigorous  solution  in  some  special  cases  [22].  The  latter  work  suggests  that, 
physically,  corner  diffraction  manifests  itself  mainly  in  the  appearance  of  “edge  wave  currents”  - 
waves  emanating  from  the  edges  and  propagating  along  the  two  edges  forming  the  corner.  We  also 
observe  such  behavior  in  our  MoM  solutions  for  the  currents. 

We  plan  to  analyze  the  problem  of  corner  diffraction  more  thoroughly,  and  we  expect  that, 
as  more  complete  calculations  of  corner  diffraction  coefficients  become  available,  they  will  improve 
the  agreement  of  the  asymptotic  theory  with  the  rigorous  results. 

Diffraction  on  a  system  of  rectangular  plates 

We  discuss  here  our  preliminary  results  for  a  system  of  two  square  plates:  the  plate  A  of  size 
10  A  x  10  A,  and  the  plate  B  of  size  8  A  x  8  A,  located  one  above  the  other,  as  shown  in  Fig.  30.  The 
plates  are  illuminated  by  a  vertically  incident  v-polarized  wave  (with  the  electric  field  along  the 
x-axis). 
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Fig.  30:  A  system  of  two  plates  (A  and  B)  illuminated  from  above  by  an  incident  plane  WF ,  so 
that  plate  B  is  in  the  shadow  of  plate  A. 


For  this  configuration,  our  present  WF  algorithm  includes,  as  scattering  mechanisms, 

(a)  reflection  on  plate  A, 

(b)  diffraction  on  plate  A, 

(c)  diffraction  on  plate  A,  reflection  on  plate  B, 

(d)  diffraction  on  plate  A,  diffraction  on  plate  B. 

We  evolve  the  WF  through  two  relatively  long  steps.  The  first  step  creates  WFs  reflected  and 
diffracted  from  the  plate  A  (Fig.  31a).  During  the  second  step  the  diffracted  wave  arrives  at  the 
plate  B  and  generates  there  secondary  reflected  and  diffracted  WFs  (Fig.  31b). 


Fig.  31:  WFs  for  the  two-plate  system  after  the  first  (a)  and  second  (b)  evolution  steps.  For  clarity, 
the  scatterer  and  the  WFs  are  cut  along  the  symmetry  plane,  and  only  one-half  of  the  configuration 
is  shown. 


The  scattered  field  is  obtained  by  integrating  the  currents  induced  on  both  plates.  The  con¬ 
tribution  to  the  cross-section  resulting  from  the  current  on  the  plate  B  includes  then  effectively 
scattering  mechanisms  (c)  and  (d)  listed  above  ((c)  expected  to  be  dominant).  Since  the  plate  B  is 
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inclined,  these  processes  are  expected  to  generate  a  certain  asymmetry  in  the  cross-section,  more 
specifically,  an  enhancement  around  6  —  30°  due  to  reflection  off  the  plate  B.  Indeed,  some  en¬ 
hancement  is  seen  in  Fig.  32,  where  we  compare  the  WF  computation  with  the  MoM  result  for  plate 
A  only  and  for  both  plates.  However,  compared  to  the  rigorous  (MoM)  result,  the  enhancement  is 
definitely  too  small. 
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Fig.  32:  The  differential  bistatic  cross  section  (in  dB,  normalized  to  A2)  for  a  single  plate  (plate 
A)  and  the  system  of  two  plates,  computed  by  using  the  Method  of  Moments  (MoM)  and  our  WF 
evolution  algorithm. 

In  Fig.  33  we  compare  current  distributions  on  the  plates,  obtained  using  the  MoM  and  our  WF 
algorithm.  Currents  on  the  plate  B  show  symmetric  and  rather  intricate  distributions;  however, 
the  WF  computation  definitely  predicts  too  small  currents  there.  That  fact  may  be  due  to  the  lack 
of  corner  diffraction  in  our  present  WF  algorithm. 

We  also  note  that  the  WF  result  shows  no  effect  of  plate  B  on  the  current  distribution  on  A 
(while  some  effect  is  seen  in  the  MoM  current  distribution).  The  reason  is  that  in  this  particular 
computation  the  WFs  reflected  and  scattered  from  plate  B  have  not  yet  reached  plate  A  (Fig.  31). 
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Fig.  33:  Distributions  of  the  absolute  values  of  the  induced  current  on  the  plates  A  and  B,  computed 
using  the  MoM  (a)  and  our  present  WF  evolution  algorithm  (b). 


5.  Algorithm  elements  requiring  further  development 

When  describing  our  WF  algorithm  we  also  mentioned  its  current  inadequacies  and  shortcomings. 
Here  we  discuss  improvements  that  can  be  introduced  in  the  further  algorithm  development. 

1.  Detection  of  diffraction  edges.  As  we  described  in  Section  2,  the  present  algorithm  identifies 
diffraction  edges  by  finding  intersections  of  rays  with  scatterer  faces.  This  procedure  is  not 
applicable  in  the  case  when  the  rays  are  exactly  tangential  to  the  (flat)  scatterer  surface  — 
a  situation  arising  in  multiple  diffraction  on  a  flat  screen  or  on  a  polyhedral  scatterer.  The 
present  algorithm  may  also  miss  edges  of  very  thin  wedges,  it  the  wedge  width  is  less  than  the 
typical  ray-ray  spacing  h. 

A  remedy  to  these  problems  is  to  preprocess  the  scatterer  geometry  and  mark  all  edges  that 
may  be  sources  of  diffraction.  Then,  in  the  WF  evolution  process  we  can  check  proximity 
of  the  rays  to  the  marked  edges,  by  effectively  testing  ray  intersections  with  cylinders  built 
around  the  edges  (the  cylinder  radius  being  comparable  to  the  minimal  ray-ray  spacing). 
This  algorithm  modification  corrects  also  the  problem  of  missing  small  end  parts  of  the  WF 

(Pig.  15)- 

2.  Detection  of  corners.  The  present  algorithms,  both  for  reflections  and  diffraction,  detects  edges 
of  scatterer  faces  by  testing  intersections  of  rays  with  faces,  and  creates  new  rays  hitting  the 
detected  edges  by  interpolating  between  existing  rays  along  the  edges  of  the  WF  mesh.  While 
this  procedure  works  well  for  sufficiently  smooth  edges,  it  may  result,  in  the  vicinity  of  sharp 
corners,  in  localization  errors  of  order  of  the  ray-spacing. 

This  algorithm  deficiency  can  also  be  corrected  by  localizing  comers  in  the  geometry  prepro¬ 
cessing  stage,  and  then  testing  ray  intersections  with  spheres  centered  at  the  corners. 
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3.  Multiple  reflections  and  diffractions  in  a  single  evolution  step.  The  present  algorithm  implements 
partially  multiple  ray  reflections  in  a  single  evolution  step,  but  does  not  allow  computation  of 
induced  currents  in  this  situation.  Neither  does  it  allow  multiple  diffraction  in  a  single  step. 
The  difficulty  in  algorithm  implementation  stems  from  the  fact  that  neighboring  rays  may 
undergo  different  numbers  of  reflections  or  diffractions,  i.e.,  have  different  evolution  histories; 
such  as  situation  precludes  straightforward  field  interpolation  required  in  evaluation  of  the 
current. 

To  correct  this  deficiency,  we  plan  to  modify  the  description  of  ray  evolution  between  the 
consecutive  WFs:  instead  of  evolving  individual  rays,  we  will  evolve  ray  tubes,  a  ray  tube 
being  a  triplets  of  rays  constituting  vertices  of  a  WF  mesh  face.  In  a  multiple  reflection  or 
diffraction  process  these  ray  tubes  will  be  split,  until  all  rays  in  a  given  tube  have  the  same 
evolution  history,  and  interpolation  becomes  possible. 

4.  Continuity  of  the  diffracted  WF  generated  in  consecutive  evolution  steps.  In  the  present  algorithm 
implementation  there  may  appear  gaps  between  the  segments  created  in  consecutive  evolution 
steps  (Fig.  16). 

This  problem  can  be  also  corrected  by  describing  ray  evolution  in  terms  of  ray  tubes  rather  than 
individual  rays.  In  this  context  a  ray  tube  provides  an  additional  connectivity  information, 
and  “memory”  preserved  between  the  evolution  steps. 

Other  algorithm  modifications  are  extensions  to  include  additional  scattering  mechanisms: 

5.  Corner  and  tip  diffraction.  These  diffraction  mechanisms  can  be  implemented  in  conjunction 
with  the  corner  detection  procedure  mentioned  in  point  2.  However,  as  we  discussed  in  Section 
2,  the  corner  and  tip  diffraction  coefficients  are  presently  available  only  for  some  special  cases, 
or  in  approximate  forms. 

6.  Gap  and  crack  diffraction.  Although  these  processes  can  be,  in  principle,  described  as  multiple 
diffraction,  it  will  be  more  convenient  and  economical  to  model  them  as  separate  diffraction 
mechanisms  characterized  by  specific  diffraction  coefficients  [23]. 

7.  Smooth-surface  diffraction.  Diffraction  on  smooth  surfaces  constitutes  a  major  extension  of 
the  present  algorithm.  We  have  developed  a  general  idea  of  the  smooth-surface  diffraction 
algorithm,  involving  splitting  of  “ray  tubes”  incident  on  the  scatterer  surface  near  the  shadow 
boundary.  Its  implementation  would  be,  however,  rather  complex. 

8.  Treatment  of  caustics.  Problems  associated  with  caustics  in  the  scattered  field  have  been  partly 
alleviated  by  utilizing  the  procedure  of  evaluating  radar  signatures  from  surface  currents. 
Caustics,  however,  may  also  appear  on  the  surface  of  the  scatterer.  In  that  case  it  might  be 
necessary,  in  the  last  diffraction  process,  to  avoid  using  the  stationary-point  approximation 
altogether,  and  describe  the  diffraction  process  in  terms  of  one  of  “incremental”  diffraction 
theory  approaches  [24, 25, 23, 26, 27]. 

In  addition,  performance  of  the  present,  preliminary,  algorithm  implementation  is  inadequate,  and 
has  to  be  improved  in  at  least  two  areas: 

9.  WF  mesh  complexification  and  simplification.  Improvement  of  this  algorithm  is  relatively 
straightforward.  It  requires  performing  operations  on  the  WF  mesh  in  a  way  as  local  as 
possible,  i.e.,  only  in  regions  where  the  mesh  has  been  modified  or  has  to  be  modified. 
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10.  Current  computation.  Efficiency  requires  that  the  present  algorithm  (Sections  3.2  and  3.3)  be 
thoroughly  revised.  We  plan  to  reverse  the  procedure  by  finding  current  points  contained  in 
specific  ray  cells,  rather  than  finding  ray  cells  containing  specified  current  points. 


6.  Assessment  of  the  Phase  I  results 

As  the  main  result  of  the  Phase  I  effort,  we  developed  and  constructed  preliminary  implementations 
of  the  following  algorithms: 

1.  Algorithm  for  free-space  WFs  evolution. 

2.  Algorithm  for  creating  WFs  according  to  the  laws  of  Geometrical  Optics  (GO),  and  for  eval¬ 
uating  fields  due  to  GO  processes. 

3.  Algorithm  for  creating  geometry  of  the  WFs  due  to  edge  diffraction  according  to  prescriptions 
of  the  Geometrical  Theory  of  Diffraction 

4.  Algorithm  for  evaluating  fields  due  to  edge  diffraction,  according  to  a  modified  version  of 
Uniform  Geometrical  Theory  of  Diffraction  (UTD). 

5.  Algorithm  for  direct  evaluation  of  scattered  fields  in  terms  of  WF  fields. 

6.  Algorithm  for  evaluating  surface  currents  based  on  the  computed  WF  field  values. 

The  algorithms  for  geometrical  WF  construction  are  applicable  in  both  frequency  and  time  do¬ 
mains.  Our  current  implementation  of  UTD  field  computation  algorithms  has  been  carried  out  in 
frequency  domain;  it  can  be,  however,  extended  to  time  domain  in  a  rather  straightforward  way 
by  incorporating  the  time-domain  UTD  being  developed  at  the  ElectroScience  Laboratory  of  Ohio 
State  University  [28, 29, 31, 32], 

The  modification  of  UTD  mentioned  in  point  4,  in  its  present  implementation,  provides  the 
correct  behavior  of  diffracted  fields  at  all  distances  (including  those  small  compared  with  the  wave¬ 
length)  for  screen-edge  problems.  By  a  correct  behavior  we  mean  here  coinciding  with  the  exact 
solution  for  the  corresponding  canonical  problem,  in  this  case  the  half-plane.  The  construction  can 
be  generalized  in  a  straightforward  way  to  the  perfectly  conducting  wedge  problem  [1]. 

As  another  comment,  relevant  to  the  point  6,  we  mention  that  computation  of  surface  currents 
may  either  yield  the  scattered  fields  (in  a  more  general  and  robust  way  than  the  direct  evaluation 
mentioned  in  point  5),  or,  alternatively,  it  can  be  used  to  construct  a  solution  Ansatz  to  be 
subsequently  used  in  solving  rigorous  asymptotic  high-frequency  integral  equations  [2, 3].  We  stress 
that  in  the  latter  application  only  the  geometrical  construction  of  WF  and  rays  is  involved,  and 
we  avoid  difficulties  associated  with  the  computation  of  field  values  (such  as  caustics) . 

The  description  of  our  WF  approach,  presented  in  this  report,  indicates  that  many  aspects  of 
the  algorithm  are  fairly  complex  both  in  formulation  and  in  implementation.  Also,  as  we  discussed 
in  the  preceding  Section,  many  algorithm  elements  require  further  development  and  substantial 
modifications.  Nevertheless,  the  results  obtained  confirm  the  two  most  important  features  of  our 
approach: 

(i)  correct  high-frequency  scaling  (number  of  rays  and  thus  computational  cost  independent  of 
frequency);  and 
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(ii)  algorithm  ability  to  treat  multiple  reflections  and  diffraction  in  a  consistent  way  (in  agreement 
with  reciprocity). 

In  view  of  these  findings,  we  believe  that  the  Phase  I  results  do  provide  a  solid  basis  for  development 
of  an  efficient  software  package  capable  of  consistent  and  accurate  modeling  of  a  wide  variety  of 
high-frequency  electromagnetic  scattering  processes.  It  would  constitute  a  significant  improvement 
compared  to  presently  available  codes  which  typically  do  not  incorporate  features  (i)  and  (ii)  just 
mentioned. 
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1 .  Introduction 


This  report  provides  a  brief  summary  of  the  technical  support  made  available  by  the 
Ohio  State  University  ElectroScience  Lab  (OSU-ESL),  to  Monopole  Research  in  their 
development  of  a  Wavefront  Evaluation  (WE)  Algorithm  for  the  efficient  ray  based 
analysis  of  electromagnetic  (EM)  radiation  and  scattering  from  electrically  large  objects 
of  practical  interest. 

Part  of  these  contributions  involved  providing  Monopole  Research  with  some  basic 
and  important  ray  field  expressions  based  on  the  uniform  geometrical  theory  of 
diffraction  (UTD)  [1-3],  and  also  providing  associated  computer  subroutines  for 
calculating  these  field  expressions  within  the  WE  algorithm.  The  UTD  solutions  which 
have  been  provided  are  in  the  frequency  domain  (FD),  as  well  as  in  the  time  domain 
(TD).  Thus,  both  FD-UTD  and  TD-UTD  type  solutions  and  their  associated  subroutines 
have  been  included  in  this  work.  Additionally,  two  new  TD-UTD  solutions  have  been 
developed  which  are  expected  to  be  particularly  useful  in  electromagnetic  compatibility 
(EMC)  and  electromagnetic  pulse  (EMP)  applications;  however,  the  computer 
subroutines  for  calculating  the  TD-UTD  fields  in  these  two  new  cases  have  not  been 
written  as  yet  in  a  user  friendly  manner.  It  is  hoped  that  the  latter,  as  well  as  additional 
FD-UTD  and  TD-UTD  solutions  of  importance  (some  of  which  have  not  been  developed 
to  date)  could  be  included  within  the  WE  framework  in  the  future  phases  of  this  study. 

The  WE  algorithm  [4],  which  is  being  developed  by  Monopole  Research,  will  allow 
one  to  track  ray  fields  more  efficiently  than  via  other  commonly  used  ray  shooting 
methods.  This  is  possible  because  in  the  WE  approach,  one  simultaneously  tracks  a  grid 
of  rays  rather  than  tracking  a  single  ray  at  a  time.  The  latter  is  possible  since  the  WE 
method  exploits  the  intimate  relation  between  a  family  of  rays  and  their  associated 
wavefront.  In  a  numerical  sense,  a  discrete  set  of  rays  can  be  employed  to  approximately 
define  a  wavefront  patch,  and  vice  versa.  Once  an  initial  wavefront  patch  is  chosen  on  the 
illuminating  or  incident  field,  one  can  then  track  this  initial  wavefront  (or  phasefront) 
patch  as  it  evolves  via  reflection  and/or  diffraction  from  the  radiating  /scattering  object, 
to  arrive  at  observer  locations  of  interest.  The  completion  of  the  EM  solution  to  the 
radiation,  scattering  and  diffraction  of  waves  from  electrically  large  complex  objects  via 
the  WE  approach  then  requires  one  to  assign  field  values  to  the  discrete  rays  defining  a 
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wavefront  patch,  as  the  latter  evolves  in  a  given  problem  configuration.  The  set  of 
wavefront  patches  defining  the  part  of  the  incident  or  excitation  wavefront  that  interact 
with  the  radiating/scattering  object  must  satisfy  certain  connectivity  requirements.  In  this 
regard,  it  is  noted  that  the  wavefront  patches  could  split  at  sharp  edges  due  to  the  process 
of  diffraction.  The  sum  total  of  the  fields  associated  with  all  the  discrete  set  of  rays  on 
each  of  the  patches  as  they  evolve  finally  provides  the  required  solution.  Of  course  one 
may  need  to  interpolate  field  values  between  the  discrete  set  of  adjacent  rays  in  order  to 
have  a  knowledge  of  the  field  at  any  point  on  the  wavefront  patch  or  patches.  Since  a 
discussion  of  the  WE  approach  may  be  found  in  [4],  the  present  report  will  focus  on 
briefly  describing  some  of  the  types  of  ray  fields  that  need  to  be  coupled  into  the  WE 
algorithm  during  the  initial  phases  of  its  development.  Other  types  of  ray  fields  may  be 
included  during  later  phases  of  this  study  to  make  the  WE  algorithm  more  versatile. 
Specifically,  the  type  of  ray  fields  considered  for  use  in  the  initial  WE  algorithm  being 
prepared  by  Monopole  Research  are  those  that  are  needed  to  account  for  some  of  the 
basic  first  order  diffraction  mechanisms  that  are  generally  present  in  large  complex 
structures  which  can  be  built  up  from  a  combination  of  smooth  convex  surfaces  and 
edges.  The  ray  field  expressions,  and  their  associated  computer  subroutines  which  have 
been  delivered  to  Monopole  Research  for  calculating  such  ray  fields  are  based  on  the 
UTD  [1-3]  as  stated  earlier.  The  UTD  extends  the  original  geometrical  theory  of 
diffraction  (GTD)  [5]  to  avoid  the  singularities  of  the  GTD  fields  at  the  ray  optical 
shadow  boundaries.  The  accuracy  and  efficiency  of  the  FD-UTD  has  been  established 
over  the  last  two  decades  via  numerous  scale  model  measurements  and  other  independent 
approaches.  On  the  other  hand,  the  TD-UTD  is  a  recent  development  [6,7].  A  TD-UTD 
solution  is  generally  obtained  by  inverting  the  corresponding  FD-UTD  solution  which  is 
assumed  to  be  available.  The  analytical  inversion  of  the  FD-UTD  into  the  TD-UTD  is 
accomplished  via  the  Analytic  Time  Transform  (ATT)  [6,9]  which  has  many  useful 
properties  that  make  it  particularly  attractive  for  inverting  ray  solutions  from  the 
frequency  to  time  domain.  The  latter  properties  of  the  ATT  are  also  discussed  in  [6,7]. 

The  motivation  for  the  TD-UTD  development  stems  from  several  factors  some  of 
which  are  enumerated  as  follows.  Exact  TD  solutions  are  available  only  for  a  handful  of 
relatively  simple  radiating  objects/configurations.  Furthermore,  it  is  more  natural  to  study 
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transient  EM  wave  phenomena  directly  in  the  time  domain  as  opposed  to  a  less  physically 
appealing  and  less  efficient  numerical  inversion  of  the  frequency  domain  solutions  into 
the  time  domain  via  the  Fast  Fourier  Transform  (FFT)  Algorithm.  Since  TD-UTD 
solutions  employ  the  same  rays  as  the  FD-UTD,  it  retains  and  therefore  exhibits  the  same 
advantages  as  the  well  established  FD-UTD;  namely,  the  TD-UTD  is  expected  to  remain 
relatively  simple  to  use  and  also  provide  a  simple  ray  picture  for  EM  radiation  and 
scattering.  Hence,  the  TD-UTD  has  the  potential,  when  developed  to  the  same  extent  as 
the  FD-UTD,  to  efficiently  solve  large,  complex  problems  for  which  exact  analytical 
solutions  may  not  be  possible.  While  the  FD-UTD  offers  localization  of  wave  effects  in 
space  via  the  use  of  rays,  the  TD-UTD  offers  localization  of  wave  effects  in  both  space 
and  time  via  space-time  rays  which  traverse  the  same  paths  as  the  FD-UTD  rays.  The 
latter  could  make  the  TD-UTD  approach  useful  in  target  identification  (ID)  application  as 
well;  in  this  regard,  one  notes  that  the  time  domain  data  needed  for  target  ID  studies  is 
expected  to  be  less  affected  by  noise  than  the  frequency  domain  data. 

The  FD-UTD  solutions  (and  associated  computer  subroutines)  which  have  been 
provided  for  integration  into  the  WE  algorithm  in  order  to  describe  some  basic  and 
important  wave  diffraction  mechanisms  are  enumerated  in  section  2.1.  All  the 
corresponding  TD-UTD  solutions  (and  their  associated  computer  subroutines)  ,  which 
have  also  been  provided  for  integration  into  a  time  domain  version  of  the  WE  analysis  of 
the  radiation  and  scattering  of  EM  waves  from  large  complex  structures,  are  enumerated 
next  in  section  2.2.  It  is  noted  that  computer  subroutines  for  all  but  two  of  the  TD-UTD 
mechanisms  are  not  available  for  integration  into  the  WE  algorithm  at  this  time,  because 
the  development  of  these  latter  two  TD-UTD  solutions  have  only  recently  been 
completed  as  discussed  in  section  3  dealing  with  new  TD-UTD  solutions.  The  computer 
subroutines  for  the  latter  two  new  TD-UTD  solutions  will  be  submitted  to  Monopole 
Research  in  future  phases  of  this  study  since  they  are  presently  not  in  a  user  friendly 
form. 

An  e+j<ot  time  convention  is  assumed  and  suppressed  for  all  the  FD-UTD  fields 
discussed  below. 
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2.  Brief  Description  of  FD-UTD  and  TD-UTD  Fields  for 
Implementation  in  the  initial  WE  Algorithm 

In  general,  a  UTD  ray  analysis  of  the  radiation/scattering  from  complex  structures 
gets  broken  down  into  a  set  of  simpler  events  such  as  incidence,  reflection  and  diffraction 
which  may  be  calculated  in  terms  of  fields  that  propagate  along  ray  paths,  in  a  sequential 
manner,  from  the  source  to  the  observer  via  these  events  on  the  complex  structure.  The 
actual  incident,  reflection  and  diffraction  events  are  often  referred  to  as  occurring  at  flash 
points  on  the  radiating/scattering  object.  To  a  first  order  of  ray  interactions,  these  flash 
points  are  points  of  reflection  and  diffraction  on  the  structure;  these  exist  in  addition  to 
direct  rays  from  the  source  to  the  observer.  It  is  assumed  in  the  present  work  that  the 
structure  is  impenetrable  and  perfectly-conducting.  Other  multiple  ray  interactions  such 
as  multiple  reflection,  reflection-diffraction  effects,  and  multiple  diffraction  effects  can 
also  be  calculated  within  the  UTD  framework;  however,  a  discussion  of  the  latter  is  not 
presented  here  for  it  is  addressed  in  [10]  for  the  FD-UTD.  Multiple  effects  (other  than 
multiple  reflection)  have  not  been  analyzed  yet  within  the  TD-UTD  framework  and  hence 
such  a  study  also  forms  a  topic  of  future  investigation  to  elevate  TD-UTD  to  the  same 
level  of  applicability  as  the  FD-UTD.  The  points  of  reflection  that  are  considered  in  this 
work  are  associated  with  such  points  when  they  occur  on  smooth  surfaces  that  may 
terminate  at  edges,  while  diffraction  points  are  allowed  to  be  on  edges  and  at  grazing 
incidence  on  smooth  convex  surfaces,  respectively. 
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2. 1  FD-UTD  for  Edge  and  Surface  Diffraction 

As  stated  above,  a  brief  summary  of  the  relevant  FD-UTD  based  ray  field  expressions 
that  are  to  be  incorporated  initially  into  the  WE  algorithm  is  given  here.  These 
expressions  allow  one  to  solve  the  radiation/scattering  of  EM  waves  by  objects  that  can 
be  built  up  from  a  contribution  of  arbitrary  smooth  convex  surfaces,  and  curved  surfaces 
containing  arbitrary  curved  edges.  Many  practical  shapes  fall  into  this  category  of  objects 
which  can  be  modeled  as  a  collection  of  curved  surfaces  and  edges  (e.g.  aircraft,  missile, 
satellite  structures,  etc.).  Other  important  effects  which  can  involve  diffraction  from  tips 
(or  comers),  coupling  between  edge  and  smooth  surface  diffraction,  double  diffraction, 
diffraction  by  non-conducting  surfaces  etc.  also  need  to  be  considered;  however,  the  latter 
can  be  incorporated  during  the  future  phases  of  the  WE  algorithm  development. 

First,  the  edge  diffraction  effect  is  summarized,  and  then  the  smooth  convex  surface 
diffraction  effect  is  treated  likewise. 

2.1.1  FD-UTD  for  Edge  Diffraction 

The  geometry  for  edge  diffraction  is  shown  in  Fig.l,  where  a  source  illuminates  an 
arbitrary  curved  wedge  which  is  assumed  here  to  be  perfectly  conducting. 


Fig.l  :  Diffraction  by  an  Arbitrary  Curved  PEC  Wedge 


In  Fig.l,  the  pattern  of  the  source/antenna  is  assumed  to  be  relatively  slowly  varying 
at  the  point  of  edge  diffraction.  A  corresponding  ray  tube  picture  illustrating  how  a  small 
incident  wavefront  patch  can  transform  into  a  edge  diffracted  wavefront  patch  is  shown 
in  Fig.2.  The  edge  diffracted  rays  lie  on  the  Keller  cone  of  diffracted  rays  as  is  well 
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known  [5,1].  The  Keller  law  of  edge  diffraction  states  that  the  half  angle  of  the  Keller 
cone  of  edge  diffraction  at  any  point  of  edge  diffraction,  Qe  exactly  equals  the  angle  that 
the  incident  ray  makes  with  the  edge  tangent  at  Qe.  Thus,  f30  =  as  shown  in  Fig.2. 


The  electric  field  at  a  point  P  anywhere  on  the  Keller  cone  at  a  distance  sd  from  Qe, 
where  ksrf»  1  and  k  =  wave  number  of  free  space  which  surrounds  the  curved  wedge,  is 
given  by 

E‘,(p)*Ei{QE)-=DA{sd)  e'jksd  (1) 

In  (1),  the  Ed(P)  and  E‘(Qe )  represent  the  edge  diffracted  electric  field  at  P,  and 
the  incident  electric  field  at  Qe,  respectively.  The  UTD  coefficient  for  edge  diffraction  is 

denoted  by  D  [1,3]  and  it  depends  on  the  frequency  of  the  incident  wave  as  well  as  the 
angles  of  incidence  and  diffraction,  and  also  the  incident  wave  polarization  as  discussed 
in  [1-3].  4/)  is  the  spatial  spread,  (or  divergence)  factor  of  the  edge  diffracted  ray 
(from  Qe  to  P)  which  indicates  how  the  energy  spreads  in  space  after  diffracting  from  the 
edge.  a{sj  )  can  also  be  obtained  from  a  conservation  of  power  in  a  diffracted  ray  tube; 
this  conservation  law  is  valid  outside  the  incident  and  reflected  ray  shadow  boundary 
transition  regions  [  ].  The  D  contains  Fresnel  transition  functions  which  not  only  keep 
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D  bounded,  but  also  provide  the  D  with  the  appropriate  behaviour  to  keep  the  total  high 
frequency  field  continuous  at  the  shadow  boundaries  where  the  geometrical  optics 
incident  and  reflected  fields  become  discontinuous  The  latter  discontinuity  results  from 
the  shadow  boundaries  for  the  incident  and  reflected  rays,  which  are  created  by  the 
presence  of  the  edge  (in  a  perfectly  conducting  or  non-penetrable  surface). 


Fig.3  :  Incident  and  Reflection  Shadow  Boundaries  ISB  and  RSB,  respectively  for  diffraction  by  a 
PEC  wedge  illuminated  by  a  plane  wave 


Fig.3  illustrates  these  shadow  boundaries  for  a  simple  wedge  configuration  when  it  is 
illuminated  by  an  incident  plane  wave.  A  boundary  layer  exists  adjacent  to  the  shadow 
boundaries  and  it  is  referred  to  as  a  transition  region  whose  size  (or  angular  extent  on 
either  side  of  the  boundaries)  depends  on  the  frequency  as  well  as  the  source  and 
observer  locations[l-3].  The  transition  region  becomes  smaller  with  increase  in 

frequency.  The  Fresnel  integrals  within  D  play  an  important  role  within  the  incident  and 

reflection  shadow  boundary  transition  regions  to  keep  D  uniformly  valid  across  these 
shadow  boundaries. 

In  (1),  the  sd  denotes  the  distance  from  Qe  to  P  (i.e.  sd  =  QEP  )  as  in  Fig.2.  If  the 

edge  is  curved  so  that  the  edge  radius  of  curvature  is  negative  (concave  rather  than 
convexly  curved  edge),  then  the  diffracted  rays  can  go  through  a  diffracted  ray  caustic 
before  arriving  to  P,  Also,  the  field  incident  on  the  edge  at  Qe  need  not  arrive  from  a 
distant  antenna  (here  shown  by  a  point  source  with  a  pattern),  but  more  generally  it  can 
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represent  an  incident  wavefront  (created  possibly  by  an  earlier  reflection  or  diffraction 
elsewhere  in  the  case  of  a  previous  interaction  before  arriving  to  Qe).  Consequently,  the 
incident  wavefront  could  arrive  at  Qe  after  passing  through  one  or  more  caustics.  For  the 
case  when  there  is  a  diffracted  ray  caustic  between  Qe  and  P,  the  A(sd )  becomes 

/  (2) 

where  n  is  the  number  of  caustics  between  Qe  and  P  with  j  =  4-1 .  Typically,  nd  =  1. 

Also,  e~Jks*  is  the  exponential  phase  delay  along  the  diffracted  ray  propagation  path  from 

Qe  to  P.  Since  the  details  on  the  derivation  of  D  in  (1)  are  available  in  [1],  no  further 

discussion  on  the  D  is  provided  here.  It  is  noted  that  the  subroutines  required  to  calculate 

D  for  edge  diffraction  have  been  made  available  to  Monopole  Research.  Finally,  it  is 
noted  that  the  magnetic  field  along  the  diffracted  ray  at  P  can  be  found  via  local  plane 

wave  conditions,  H  (P)  ~Y0sd  xEd 

where  Yo  -  Zo'1  and  Zo~  intrinsic  free  space  wave  impedance. 

2.1.2  FD-UTD  for  Slope  Edge  Diffraction 


Fig.  4  :  Slope  Diffraction  by  an  Arbitrary  Curved  PEC  Wedge 


When  the  incident  field  is  rapidly  varying  at  any  point  of  edge  diffraction,  Qe  as 
shown  in  Fig.  4,  then  the  result  in  (1)  for  Ed(p)  must  be  augmented  by  a  slope 
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diffracted  field  contribution  Esd(p).  The  Esd(p)  is  required  to  make  the  slope  of  the 
total  field  continuous  across  incident  and  reflection  shadow  boundaries, 
because  Ed(P)  is  valid  only  for  non-rapidly  varying  incident  fields  at  Qe-  One  notes  that 
(1)  can  be  expressed  in  matrix  form  as 


-D. 


0 


0  -D, 


E‘ , 

Po 

E‘ 


Mf) 


-Jks 


(3) 


where  the  incident  and  diffracted  fields  are  expressed  in  terms  of  unit  vectors  fixed  in  the 
incident  and  diffracted  rays,  and  in  particular  in  the  edge  fixed  planes  of  incidence  and 

diffraction  as  defined  in  [1-3].  One  employs  these  special  set  of  unit  vectors  so  that  D 
then  reduces  to  a  two  term  dyad  [1].  Thus, 


Ed{P)=  PoEi  +  <fe 


E'(P)=0,E‘t.+}  4  ; 


A  t 


*  A  _  1 


(4) 

(5) 


D  =  -fi9  fi0  D.-M  A,  (6> 

One  notes  that  the  are  perpendicular  to  each  other  and  to  sd ,  while 

are  perpendicular  to  each  other  and  to  the  incident  ray  direction,  s' . 

The  slope  diffracted  field  Esd{P)  for  a  perfectly  conducting  wedge  can  be  expressed 
more  conveniently  in  matrix  form  as 


rpsd 

Dsi  —E*.  DSI  9  Er. 
dri  P*  dnr  P' 

EPo 

Ef_ 

Dsi  9 .  E\  D?  9  Er. 
dn'  *  "  dnr  *r  _ 

(7) 


Again,  the  unit  vectors  fiQ  ,  and  ^  are  as  defined  in  [1],  whereas,  and  <j>'r 

are  the  corresponding  unit  vectors  fixed  on  the  reflected  ray  and  in  the  edge  fixed  plane 
of  reflection  as  defined  in  [1].  Also,  the  slope  diffraction  coefficients  If,  Efr  and  Df  are 
defined  in  [2];  they  contain  slope  transition  functions  that  are  of  the  type 

2  jx(\  -  F{x))  (8) 
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where  the  ordinary  transition  function  F(x)  occurs  directly  in  the  ordinary  UTD  edge 

diffraction  coefficient  D  (or  in  Ds  and  A,,  respectively ).  Here  F(x)  contains  the  Fresnel 
integral  alluded  to  earlier;  in  particular, 

oo 

F(x)=2jyfx  ejx  J  e~Jfl  dr  (9) 

Within  the  transition  region  x  is  small  and  close  to  the  shadow  boundary  so 
F(x)  « *Jnjx ,  whereas  outside  the  transition  region  x  is  large  (usually  x>3)  so  that 

F(x)  =  1 .  More  details  on  this  slope  diffracted  field  E sd  may  be  found  in  [2]  and  also  in 
Appendix  1 .  Also,  the  computer  subroutines  for  calculating  the  slope  transition  function 
of  (8)  and  hence  the  UTD  slope  diffraction  coefficient  have  been  made  available  to 
Monopole  Research.  Finally,  one  notes  that  the  slope  diffracted  field  allows  for  the 
incident  and  reflected  fields  to  vary  rapidly  along  the  directions  normal  to  the  edge  fixed 
planes  of  incidence  and  reflection,  respectively;  i.e.  along  ft'  and  n  as  is  evident  from 
the  normal  derivatives  given  by  h'  •  V  and  hr  ■  V  in  (7). 

2.1.3  FD-UTD  for  Surface  Diffraction 


Fig.  5(a)  :  Diffraction  at  a  Smooth  Convex  Surface 


When  an  incident  ray  grazes  a  smooth  perfectly  conducting  convex  surface  at  Qi,  it 
launches  a  surface  ray,  which  propagates  along  a  geodesic  path  into  the  shadow  region  of 
the  convex  object  as  shown  in  Fig.5  (a)  and  Fig.  5(b). 
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/ 

Fig.  5(b) :  Unit  Vectors  fixed  in  the  Surface  Diffracted  Ray 

Furthermore,  the  field  along  the  ray  geodesic  path  is  exponentially  attenuated  because 
energy  is  shed  from  the  surface  ray  along  its  forward  tangent.  This  energy,  which  is  shed 
from  the  surface  ray,  travels  a  surface  diffracted  ray  path  to  arrive  at  an  observer  at  Ps  in 
the  shadow  region  of  the  smooth  convex  object  as  shown  in  Figs.  5(a)  and  5(b).  The 
surface  diffracted  ray  arrives  at  P  from  the  point  of  surface  diffraction  on  geodesic 
surface  ray  at  Q2.  The  path  from  point  of  launching  of  the  surface  ray  at  the  point  of 
grazing  incidence  Qi  to  the  point  Q2  is  the  surface  ray  path  which  is  a  surface  geodesic.  A 
surface  ray  strip  and  its  associated  surface  diffracted  ray  tube  are  depicted  in  Fig.  6  in  the 
general  case  when  an  astigmatic  ray  tube  is  incident  at  grazing  on  a  smooth  convex 
surface.  If  the  smooth  convex  body  is  a  closed,  as  it  usually  is  in  practice,  the  surface  ray 
can  also  traverse  around  the  body;  hence,  it  can  produce  surface  diffracted  rays  which  can 
also  reach  an  observer  in  the  lit  region  at  Pi  where  the  incident  ray  field  is  directly 
visible.  In  the  lit  region,  the  total  field  at  Pi  consists  of  the  incident  and  reflected 
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geometrical  optics  field  as  well  as  surface  diffracted  ray  fields  for  closed  bodies.  In  the 
shadow  region,  the  total  field  at  Ps  consists  of  only  the  surface  diffracted  ray  fields. 


Fig.  6  :  Surface-diffracted  ray  tube.  Radius  of  the  incident  wavefront  where 
it  grazes  the  surface  at  Qi  is  p1 


In  the  smooth  convexly  curved  body,  the  incident  and  reflected  ray  directions  merge 
along  grazing  incidence  to  define  the  surface  shadow  boundary  (SSB).  The  SSB  is  an 
extension  of  the  incident  ray  past  grazing.  A  boundary  layer  exists  adjacent  to  SSB;  this 
layer  is  referred  to  as  the  SSB  transition  region.  On  the  lit  side  which  is  outside  this 
transition  region,  the  field  consists  of  ordinary  geometrical  optics  incident  and  reflected 
ray  fields  for  a  given  incident  ray  field  illumination.  In  addition,  one  or  more  surface 
diffracted  ray  field  contributions  can  also  be  present  in  the  lit  region  as  mentioned 
previously.  The  latter  occur  if  the  surface  rays,  which  continuously  shed  surface 
diffracted  rays,  encircle  the  body  one  or  more  times.  Effects  of  multiple  surface  ray 
encirclements  may  be  ignored  if  the  closed  convex  body  is  electrically  large,  because  the 
surface  ray  field  in  this  case  becomes  exponentially  weak  (due  to  continental  shedding  of 
energy)  to  the  point  where  it  can  be  neglected.  On  the  other  hand,  the  ordinary 
geometrical  optics  reflected  field  vanishes  along  the  SSB;  however,  there  is  a  diffraction 
effect  present  which  modifies  the  reflected  field  at  and  near  the  SSB  (i.e.  on  the  lit  side  of 
the  SSB).  Hence,  the  UTD  electric  field  at  a  point  Pi  (see  Fig.5  (a))  on  the  lit  side  of  the 
SSB  is 

E(PLh  E‘(Pl)+  E‘r(PL)+E*{PL)  (10) 
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where  Ed(PL )  results  from  single  or  multiple  encirclements,  and  can  be  ignored  if  the 
closed  convex  body  is  large  as  mentioned  before.  The  Egr(PL)->  Er(PL)  in  the  lit  region 
outside  the  SSB  transition  zone,  and  Er(PL )  is  the  conventional  geometrical  optics 
reflected  field.  Hence  E gr  (PL )  is  the  generalized  reflected  field  which  contains 
diffraction  separate  from  that  in  E d {PL )  due  to  encirclements.  One  notes  that  Egr(PL)  is 

far  more  significant  than  E  d  (PL )  in  the  lit  zone. 

At  a  point  Ps  on  the  shadow  side  of  the  SSB,  the  total  field  is  simply  the  surface 
diffracted  field  Ed(Ps),  i.e., 

E(P,)*Ed(P,)  (11) 

Additional  contributions  to  Ed{P])  from  encircling  rays  can  also  be  present  at  Ps  ; 
however;  their  contribution  is  negligible  for  large  convex  bodies. 

The  fields  E^Pj,  Egr(PL)  and  Ed(Ps)  indicated  in  (10)  and  (1 1)  are  given  in  detail 

in  [11,12].  The  form  of  Ed(PL)  is  the  same  as  Ed(Ps).  The  general  forms  of 
Egr{PL) and  Ed(Ps)  are  [11,12] : 


S'fo)- £'(&)•*  Ar{sr)  e~Jks'  (12) 

!'(/>>  £'(&)•  5  Ad(sd)  e**  (13) 


Outside  the  transition  on  the  lit  side  region,  R-*  R,  where  R  is  the  dyadic 
geometrical  optics  reflection  coefficient.  Likewise,  outside  the  transition  region  on  the 

shadow  side,  D  reduces  to  Keller’s  field  description  in  terms  of  GTD  surface  diffraction 
coefficients  at  Qi  and  Q2,  respectively  together  with  a  complex  propagation  constant 

along  the  geodesic  path  from  Qi  to  Qj  [13].  Within  the  transition  regions,  both,  £  and  D 
in  the  UTD  solution  contain  an  F  type  transition  function  as  in  (9)  and  a  Psh  ,Fock  type 
transition  function,  which  keep  the  total  field  bounded  and  continuous  across  the  SSB, 


respectively.  From  [11,12],  one  can  write  D  as 


D  = 


+  D,n,n2 


(14) 
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where  )  are  the  normal  and  binormal  unit  vectors  to  the  surface  ray  geodesic  at  the 

launch  point,  Qi,  while  («2,£2)  are  the  corresponding  unit  vectors  at  the  point  of 

shedding,  Q2  as  in  Fig.  5(b).  The  unit  vectors  )  are  along  the  forward  tangents  to 

the  surface  ray  at  at  Qi  and  Q2,  respectively.  One  notes  that  6,  =  tt  x  nx .  Also,  t  in  (14)  is 

222 

the  geodesic  surface  ray  arc  length  from  Qi  to  Q2.  Furthermore,  ,4,.(sr)  and  Ad(sd)  are 
the  spatial  spreading  (or  divergence)  factors  associated  with  the  reflected  ray  (from  Qr  to 
Pi )  and  with  the  surface  diffracted  ray  (from  Q2  to  Ps).  The  di](Q]  2)  in  (14)  represents 

the  width  of  the  surface  ray  strip  at  Qu  as  in  Fig.  6  and  hence  the  square  root  in  (14) 
represents  the  spreading  of  energy  in  a  surface  ray  strip  along  a  geodesic  ray  .  path. 

Finally,  sr  and  s<j  are  the  reflected  and  surface  diffracted  ray  distances  ( QRPL  )and 


( \Qi I )  respectively. 

It  is  noted  that  the  UTD  dyadic  transition  functions  in  (12)  and  (13)  can  be  computed 
via  the  subroutines  that  have  been  made  available  to  Monopole  Research. 


2.1.4  FD-UTD  for  the  Radiation  from  a  Source  on  a  Smooth  Convex 
Surface 

When  a  source/antenna  is  placed  on  a  smooth  convex  surface,  which  is  perfectly 
conducting,  then  the  near  and  far  field  radiated  by  this  source  at  Q[  can  be  shown  to  be 
associated  with  a  direct  ray  to  a  point  Pi  in  the  lit  region  (where  the  source  is  directly 
visible),  and  a  surface  diffracted  ray  to  a  point  Ps  in  the  shadow  region.  These  rays  are 
depicted  in  Figs.  7(a)  and  7(b),  respectively,  for  a  non  closed  smooth  convex  body.  If  the 
convex  body  is  closed,  as  is  usually  the  case  in  practice,  then  additional  oppositely 
encircling  surface  rays  can  also  contribute  via  surface  diffraction  to  the  point  Pi.  These 
additional  contributions  can  also  be  present  at  Ps.  Here,  these  additional  contributions  are 
assumed  to  be  negligible  as  is  true  for  large  convex  bodies. 
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7(a)  7(b) 

Fig.  7  :  Radiation  by  a  Source  on  a  Smooth  Convex  Surface 


As  shown  in  [14],  the  electric  fields  dEme (PL )  and  dEmje(Ps)  radiated  by  an 
infinitesimal  magnetic  (m)  or  electric  (e)  current  moment  dpme(Q')  at  Qi  can  be 


expressed  symbolically  as 


dE.,.(PL>dp.,{Q')  T'.,A^d) 


(15) 


and 

e-*  (16) 

where  A(s)  (  =  1/s)  and  Ad(sd)  constitute  the  spreading  of  the  rays  along  the  direct  path 


along  5  (=  \q'Pl  ),  and  the  surface  diffracted  path  sd  (=  | QPS ),  respectively.  One  notes 


that  the  ray  path  from  at  to  Q  (before  it  sheds  from  Q  to  Ps )  is  a  surface  geodesic. 
Also,  the  surface  ray  divergence  factor  from  to  Q  is  given  by  the  square  root  term  in 
(16);  the  dy/a  and  <f  77(g)  in  (16)  are  shown  in  Fig.  8.  The  dyadic  radiation  and  diffraction 


functions  Tm,e  and  Tm,e ,  respectively,  are  given  in  terms  of  UTD  transition  functions 
involving  the  radiation  Fock  functions  g  and  g  [14,15].  The  unit  vectors  (?' ,b' )  at  (/ 
and  (/  ,n,b)  at  Q,  respectively,  are  fixed  in  the  surface  diffracted  ray  as  in  Fig.  7(a);  the 
Tm.e  then  shows  how  a  given  orientation  of  dpme  at  Q  (i.e.  t' ,n'  or  b')  produces  the 


„  _  —r.  A  ~ 

polarization  type  ( n ,  b  )  for  dE (Ps )  at  Ps.  Likewise,  T m,e  indicates  how  the  (t'nb'(  and  h ' ) 

* 

components  of  dpme(Q')  at  Q[  produce  the  (n,b)  polarized  fields  dE(PL)  at  Pi.  Thus, 

Tm,e  and  Tm,e  are  the  UTD  transition  or  transfer  functions  (containing  the  Fock  functions 
g  ,  g )  that  link  the  source  to  the  fields  produced  by  the  source;  they  also  keep  the  total 
field  continuous  as  one  goes  from  Pi  on  the  lit  side,  to  Ps  on  the  shadow  side  of  the 
surface  shadow  boundary.  These  UTD  dyadic  transition  functions  in  (15)  and  (16)  can  be 
computed  via  the  subroutines  made  available  to  Monopole  Research. 


Fig.  8  :  Surface  Diffracted  Ray  Tube  excited  by  a  Source  at  on  a  Smooth  Convex  PEC  Surface 

The  response  Eme(PL)  and  Eme(Ps )  due  to  any  physical,  finitely  distributed 

equivalent  magnetic  (m)  or  electric  (e)  sources  can  be  obtained  by  integrating  dEm  e  (PL ) 

and  dEm  e (Ps )  over  these  sources  [14].  Although,  the  results  in  (15)  and  (16)  provide  the 

fields  radiated  by  an  infinitesimal  source  on  a  smooth  convex  surface,  they  can  also,  by 
the  reciprocity  theorem,  provide  the  surface  charge  density  (e  case)  or  the  surface  current 
density  (m  case)  at  which  is  induced  by  a  wave  incident  at  (/  from  a  distant  source  at 
Pi  or  Ps,  respectively.  The  latter  is  of  interest  in  electromagnetic  pulse  (EMP)  and 
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scattering  problems,  while  the  solutions  in  (15)  and  (16)  given  for  the  direct  problem  are 
useful  for  predicting  the  radiation  by  conformal  antennas  and  arrays  on  smooth  convex 
surfaces. 

2.1.5  FD-UTD  for  the  Mutual  Coupling  between  Antennas  on  a 
Smooth  Convex  Surface 


geodesic  surface 


•  Surface-ray  strip  (or  tube) 

Fig.  9(a) :  Surface  ray  from  Q  to  Q 


POINT] 


Fig.9(b)  :  Surface-ray  strip  (or  tube) 
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The  surface  ray  field  at  any  point  Q  on  a  smooth,  perfectly  conducting  convex  surface 
which  is  produced  by  a  source  at  (/  on  the  same  surface  must  again  follow  a  geodesic 
path  connecting  to  Q.  The  ray  geometry  for  this  problem  is  sketched  in  Figs  9(a)  and 
9(b). 

The  electric  and  magnetic  surface  fields  dEm  e  and  dHm  e  at  Q  may  be  expressed  via 
the  development  in  [16]  as  follows: 

<*., fe)  -  (S')- «'*'  0T> 

and 

dH,..(Q)*dP',.(3') (18) 

The  dy/0  and  drj{Q),  both  of  which  appear  in  (17)  and  (18),  respectively,  have  the 
same  meaning  as  in  (16);  these  quantities  are  also  depicted  in  Fig.  8(b).  The  surface  ray 

=e  =//  A  a 

field  dyadics  r«,m  and  Te,m  indicate  the  type  of  surface  field  polarization  (t  ,n  or  b ), 

A  A 

which  is  produced  at  Q  by  any  particular  component  (/',«'  or  b')  of  the  source 

—e  =h 

dpm  e\Q')  at  Q.  A  detailed  description  of  T e<„,  and  Te,m  are  available  in  [16].  The 
A(s)  =  -L  in  (17)  and  (18),  where  s  is  the  geodesic  arc  length  from  (/  to  Q  as  sketched 
in  Fig  9(a). 

The  UTD  dyadic  transfer  or  transition  functions  in  (17)  and  (18)  contain  the  surface 
Fock  functions  U  and  V,  respectively  [16,15];  these  UTD  transition  functions  can  be 
computed  via  the  subroutines  made  available  to  Monopole  Research. 
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2.1.6  An  Application  of  FD-UTD  To  Airborne  Antenna  Pattern 
Prediction 


As  indicated  previously,  the  FD-UTD  has  been  utilized  extensively  to  analyze  the 
radiation  from  complex  shapes.  The  FD-UTD  results  have  been  tested  numerous  times 
with  scale  model  measurements  or  by  other  independent  numerical  approaches  where 
possible.  An  example  analyzed  here  via  the  FD-UTD  is  indicated  below  to  illustrate  the 
versatility  of  this  ray  approach.  Specifically,  a  complex  cavity  backed  slot-blade  antenna, 
which  is  analyzed  via  the  moment  method  (MM)  solution  of  the  governing  integral 
equation  is  placed  on  F-16  fighter  jet  in  this  example.  The  pattern  of  this  antenna  on  F-16 
is  analyzed  via  the  Ohio  State  Univ.  (OSU)  NEW AIR  code  which  is  based  on  the  FD- 
UTD  [17].  The  antenna  geometry,  as  well  as  the  NEW  AIR  FD-UTD  model  of  F-16,  and 
the  patterns  calculated  from  that  code  are  illustrated  in  Figs.  10,  11  and  12,  respectively. 
It  is  important  to  note  that  the  FD-UTD  patterns  in  Fig.  12  have  been  computed  based 
primarily  on  the  FD-UTD  formulas  given  in  item  2.1.1  and  item  2.1.4  above. 
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It  would  be  of  interest,  for  the  reasons  mentioned  in  section  1,  to  develop  a 
corresponding  TD-UTD  which  would  allow  one  to  analyze  the  radiation  /scattering  of 
transient  EM  waves  from  complex  objects  to  the  same  level  as  the  FD-UTD.  The  recent 
TD-UTD  developments  are  described  next  in  sections  2.2  and  2.3. 
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Fig.  12  :  FD-UTD  Radiation  Patterns  of  the  Antenna  in  Fig.  10  when  placed  on  the  aircraft  of  Fig.  1 1 

2.2  TD-UTD  for  Edge  and  Surface  Diffraction 

The  motivation  for  the  development  of  a  TD-UTD  to  efficiently  analyze  the  transient 
radiation  /scattering  from  complex  objects  which  are  excited  by  short  pulses  has  been 
discussed  earlier  in  section  1.  The  TD-UTD  employs  the  same  rays  as  the  FD-UTD  and 
hence  retains  all  of  the  advantages  of  the  FD-UTD  since  one  constructs  a  TD-UTD 
solution  from  an  analytic  inversion  into  the  time  domain  of  the  corresponding  FD-UTD. 
The  TD-UTD  solutions  which  have  been  developed  recently,  but  prior  to  the  start  of  the 
present  contract  with  monopole  Research,  are  summarized  below.  The  particular 
transform  used  in  the  TD-UTD  development  which  converts  the  FD-UTD  into  the  time 
domain  is  the  Analytic  Time  Transform  (ATT).  The  ATT  is  utilized  because  it  has 
several  advantages  for  this  particular  application.  Specifically,  the  ATT  is  convenient  for 
inverting  the  FD-UTD  fields  which  may  have  traversed  through  ray  caustics,  and  which 
may  exhibit  a  general  elliptic  polarization. 

The  ATT  is  defined  for  positive  frequencies  (i.e.,  for  ©  >  0)  as  follows: 

/(t)  =  — ;  Im t>a  ,  (19) 

n  o 

where  /(f)  is  the  ATT  of  the  frequency  domain  function  F(co).  fit)  is  analytic  for 
Im  t>a  with  \F{coJ  *  ceao>  as  a  ->  oo .  The  f(t)  is  obtained  from  /(/)  via 
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f(t)=Ref(t)  ;  lmf  =  0  , 

Note  that  f(t)  is  the  inverse  Fourier  transform  of  F(co).  Also, 


f{t)=f(t)+jHf(t)  ; 

in  which  Hf(t)  is  the  Hilbert  transform  of f(t). 
The  analytic  delta  function  is  defined  as 


Imf  =  0 


<5(r)+  pv— 
nt 


,lmt  >  0 


,Imr  =  0 


where  pv  signifies  the  Cauchy  principal  value  when  integrating  over  the  function.  Hence 

£(f)  is  a  distribution  for  Im/  =  0 ;  it  may  be  used  as  an  input  or  excitation  with  a  time 
impulsive  behaviour  for  the  transient  TD-UTD  solution  that  is  to  be  constructed.  The  TD- 

UTD  response  to  a  S(t)  type  excitation  is  generally  easier  to  obtain;  it  then  allows  one  to 
construct  the  response  to  any  realistic  finite  bandwidth  excitation  via  an  efficient 
convolution  procedure.  For  example,  if  a  realistic  excitation  pulse  E'0{o))  can  be 
represented  in  the  frequency  (co)  domain  by  a  sum  of  handful  of  exponential  functions  as 


,o)  >  0 


&t  +  jan 


,  Imf  >  -a„ 


Now,  if  ei(t)  denotes  the  TD-UTD  response  to  an  analytic  delta  (impulse)  function 

+ 

i.e.  to  a  pS{t)  time  dependence,  then  the  TD-UTD  response,  e{t)  is  given  via 
convolution  as 

Sw-4£rr t  *  •'«  (25) 
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£  >  0 
Im/  >  a 


(26) 


N  x+je  + 


J-(''Vz7 hr*1 

Z;r  «=l  -oo+ y*  f  1  Jan 


_T 

Since  e/(/')  may  be  assumed  to  be  analytic  for  Im  t'  >  a  and  the  poles  of  —  ” 


reside  in  the  same  half  space,  it  follows  that 

^)=YjAneiit  +  jan)  (27) 

«=1 

+ 

As  long  as  Re«„  >0  for  all  n,  the  e{t)  will  be  analytic  on  the  real  time  axis 
(Imt  =  0)  so  one  can  obtain  the  desired  TD-UTD  response  e(t)  in  a  straight  forward 
manner  via 


e(t)=Ree(/)  ;  lmf  =  0  .  (28) 

Consider  a  typical  waveform  shown  in  Fig.  13  whose  frequency  (a>)  domain 
behaviour  is  shown  in  Fig.  14.  The  behaviour  in  Fig.  14  can  be  synthesized  with  only 
three  terms  of  the  type  in  (23);  consequently,  the  response  in  (27)  and  hence  that  in  (28) 
can  be  found  by  summing  (24)  over  just  three  terms.  The  attractive  property  that  allows 
one  to  go  from  (26)  to  a  simple  result  for  the  convolution  in  (27)  is  another  reason  to 
employ  the  ATT  rather  than  the  conventional  Fourier  or  Laplace  transforms  in  the  TD- 
UTD  development. 

Some  TD-UTD  results  available  from  the  recent  past  are  reviewed  next. 


Fig.  13:  The  waveform  in  Time  Domain  Fig.  14  :  Frequency  Domain  Behaviour 
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2.2.1  TD-UTD  for  Edge  Diffraction 


The  geometry  for  edge  diffraction  is  shown  in  Fig.l  as  before  for  the  FD-UTD  case. 
The  TD-UTD  field  for  diffraction  by  an  arbitrary  curved  conducting  wedge  can  be  found 

by  taking  the  ATT  of  (1)  for  a  general  astigmatic  wave  illumination.  The  case  of  the 

+ 

astigmatic  wave  illumination  in  the  FD-UTD  case  leads  to  a  time  impulsive  (S) 

astigmatic  wave  illumination  for  the  corresponding  TD-UTD  case.  Hence,  the  impulsive 

+ 

response  ef  {t)  for  the  wedge  diffraction  has  been  obtained  from  the  ATT  of  (1)  for  this 

situation  as  [6,7] 

+  +_ 

e1(>)*‘K(ac)j(rMs‘‘)  (29) 

where 


s  s 
c  c 


rd  =t - 


(30) 


in  which  s,  and  Sd  are  the  distances  from  O  to  Qe  (i.e  s'  =  \OQE\  )  and  Qe  to  P  (i.e 


sd  = 


\QeP\  38  before).  Also,  c  =  speed  of  the  EM  wave  in  free  space  which  surrounds  the 
wedge.  Since  the  details  of  the  development  leading  to  (29)  are  available  in  [  ,  ],  no 

further  details  will  be  provided  here,  except  to  state  that  the  d(rd)  of  the  TD-UTD  turns 


out  to  be  actually  simpler  than  the  corresponding  D  of  (1)  for  the  FD-UTD.  Subroutines 
for  computing  (29)  have  been  made  available  to  Monopole  Research. 


2.2.2  TD-UTD  for  Slope  Edge  Diffraction 

For  the  case  of  an  incident  wave  which  exhibits  a  rapid  spatial  variation  at  the  point 
of  edge  diffraction,  a  slope  diffracted  field  contribution  had  to  be  added  to  (1)  as 
indicated  in  2.1.2.  The  details  of  the  TD-UTD  development  from  the  corresponding  FD- 
UTD  given  in  (7)  are  also  available  in  [6],  However,  a  paper  on  this  subject  has  been 
written  on  the  present  contract  [18]  and  it  is  also  included  as  Appendix  1  of  this  report  for 
the  sake  of  completeness.  Thus,  no  further  details  about  a  TD-UTD  for  slope  edge 
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diffraction  are  provided  here.  The  computer  subroutines  for  calculating  TD-UTD  slope 
diffracted  field  have  been  made  available  to  Monopole  Research. 

2.2.3  TD-UTD  Surface  Diffraction 

The  geometry  of  this  surface  diffraction  problem  is  found  in  Fig.  5  for  the 

+ 

corresponding  FD-UTD  case.  The  TD-UTD  solution  for  a  time  impulsive  (<5)  astigmatic 
wave  excitation  can  be  found  by  taking  the  ATT  of  (10)-(13)  for  the  corresponding  FD- 
UTD  case  of  an  astigmatic  ray  optical  field  illumination  of  the  smooth  convex  conducting 
surface  as  described  in  [6],  A  paper  describing  this  TD-UTD  development  is  in 
preparation  as  a  part  of  the  present  work  for  Monopole  Research  [19].  The  desired  TD- 

UTD  solution  is  given  as 
+  + 

,  lit  side  of  SSB 

(31) 

+ 

ei(t)  ,  shadow  side  of  SSB 

where 

=  (32) 

+  _  jL 

(33) 

and 

+ 

e(r)=4  (34) 

+  + 

The  R  and  D  are  the  ATT  of  the  transition  functions  involving  the  F  and  Ps  h  type 
integrals  in  (10)-(13). 

As  mentioned  above,  the  details  concerning  the  application  of  the  ATT  leading  to 
(31)-(34)  are  available  in  [6],  and  will  soon  appear  in  [19].  Computer  subroutines  for 
calculating  (31)-(33)  have  been  made  available  to  Monopole  Research. 
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3.  Summary  of  Two  New  TD-UTD  Solutions 

The  TD-UTD  solutions  for  the  radiation  from  antennas  on  a  smooth  convex*  surface, 
and  the  mutual  coupling  between  a  pair  of  antennas  on  a  smooth  convex  surface, 
respectively,  have  been  completed  under  the  present  contract.  These  new  TD-UTD 
solutions  are  developed  via  the  application  of  the  ATT  to  the  general  expressions 
indicated  in  sections  2.1.4,  and  2.1.5,  respectively,  for  the  corresponding  FD-UTD 
solutions;  these  solutions  are  summarized  below  together  with  the  some  numerical  results 
illustrating  their  accuracy. 


3.1  TD-UTD  for  the  Radiation  from  a  Source  on  a  Smooth 
Convex  Surface 


The  geometry  of  the  problem  in  question  is  the  same  as  that  in  Fig.  7(a).  The  FD- 
UTD  for  this  problem  is  given  separately  for  the  lit  and  shadow  regions  as  in  (15)  and 
(16),  respectively. 

The  ATT  of  (15)  or  the  more  detailed  version  in  [12]  yields  the  TD-UTD  solution  for 
the  lit  side  as: 


d Em,e{t )  -  dpm  e 


=i.t 
•  T  m,e 


Mil 


V  cjs 


(35) 


with 


=i,t , , 


-1 


Am 


*  + 


b[h  A+  t[b  B+  b'eb  C+  t'th  D 


(36) 


-Zn 


Am 


h'enM+i'(fiN 


(37) 


where 

A  =  He+T02T  cos0t  (38) 

B  =  Si-Tq  T cos2 (39) 


C  =  T0  f  (4°) 

z)  =  r0f  cos  #(.  (41) 
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M  -  sin  #, 


H(  +  TqT  cos#, 


N  =  sin#i7’0  T 


with 


cos#, 
"  [i+r02cos#,' 


(42) 

(43) 

(44) 


and 


H,(S„t)=Arr[H'(a,,$,)] 

S'(3„t,M,)  =  Arr[s'(a>,Z„m,)] 
The  parameters  in  the  lit  region  are  given  by 


M(.) 


(45) 

(46) 

(47) 


[l  +  r02cos2#,f 

Z,  =-M/(*)cos#,  (48) 

One  can  also  take  the  ATT  of  (16)  the  more  detailed  version  in  [12]  to  obtain  the  TD- 
UTD  solution  for  the  shadow  side  as: 

=/  ( 


dEm,e{t)=dp°m  e  •r* 


S,  S 
k  C  C  J\S 


*(pc+S) 


(49) 


with 


T  -  _1 


4m: 


b'hT ;  (Q')H+  t'bT2  (Q')s+  b'bT3  {Q')S+  t'nT \  (Q')H 


dy o 

1/6 


(50) 


and 


-Z„ 

=  0 


4m: 


I  dy 0 

1/6 


(51) 


where 


H(Z,t)  =  ATT[H(co,Z)] 
S{Z,tM)=ATT[s(a),Z,m)] 


(52) 

(53) 


and 
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(54) 


M(.)  = 


,(•) 


\'/3 


2  c 


(55) 

The  To  in  the  above  equations  is  referred  to  as  the  torsion  factor  [12]  which 
characterizes  torsional  surface  rays. 

Algorithms  have  been  developed  to  calculate  the  ATT  of  the  relevant  radiation  Fock 
functions  g  and  g  which  appear  in  the  above  TD-UTD  solutions.  These  subroutines  are 
presently  not  written  in  a  user  friendly  form.  It  is  planned  to  make  these  subroutines  user 
friendly  in  the  future  phases  of  this  study.  The  details  of  this  new  TD-UTD  solution  will 
be  submitted  for  publication  soon  [20].  It  is  noted  that  the  above  TD-UTD  solutions  have 

been  obtained  for  an  excitation  J  or  M  which  is  d p e  m5(r  -r')u(t  - 1')  where  U  is  a 

step  function,  and  J  or  M  are  the  source  current  densities  at  f/.  The  fields  radiated  by  a 
time  impulsive  source  at  Q[  may  be  found  by  differentiating  the  step  response  which  is 
more  easily  obtained  first. 

Some  numerical  results  are  illustrated  in  Figs.  15-20  for  a  source  at  £/  on  a 
circular  cylinder,  where  the  transient  pulse  excitation  at  Q  is  the  same  as  in  Fig.  13;  The 
results  for  this  pulsed  excitation  are  found  via  an  ATT  based,  essentially  closed  form, 
convolution  as  discussed  at  the  beginning  of  section  2.2.  The  circular  cylinder  geometry 
is  chosen  because  an  exact  frequency  domain  eigenfunction  solution  can  be  constructed 
for  this  case  and  transformed  into  the  time  domain  for  comparison.  The  results  in  the  far 

A 

zone  at  a  given  angle  if)  are  shown  for  both  z  (axial)  directed  as  well  as  <j> 

(circumferentially  directed)  magnetic  current  sources  at  (x  —  a,  y  —  0t  z  =  0)  where  a  is 
the  radius  of  the  cylinder.  Here,  a  —  l  meter.  The  agreement  between  the  TD-UTD  and 
the  reference  solution  (referred  here  simply  as  eigen  solution)  is  very  good. 
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Fig.  15:  TD-UTD  Radiation  from  a  Source  on  a  Circular  Cylinder,  which  is  evaluated  in  the  far  zone  in  the 

direction  <j>. 


Fig.  16:  TD-UTD  Radiation  from  a  Source  on  a  Circular  Cylinder,  which  is  evaluated  in  the  far  zone  in  the 

direction  <|). 


Fig.  17:  TD-UTD  Radiation  from  a  Source  on  a  Circular  Cylinder,  which  is  evaluated  in  the  far  zone  in  the 

direction  <|>. 


E=iEM=  35') 


direction 


Fig.  19:  TD-UTD  Radiation  from  a  Source  on  a  Circular  Cylinder,  which  is  evaluated  in  the  far  zone  in  the 

direction  (j>. 
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Um  =  ATT 


(60) 


Vm 


(61) 


The  convolution  with  s{t  -  p/J  in  the  above  formulas  will  simply  cause  a  time 

delay  due  to  surface  ray  propagation  along  the  geodesic  path  length  s  from  ^ t0  Q • 

Algorithms  have  been  developed  to  calculate  the  ATT  of  the  relevant  surface 


+  + 

Fock  functions  U  and  V  which  appear  in  the  above  TD-UTD  solution;  however,  they 
are  not  presently  available  in  a  user  friendly  form  but  they  will  be  made  available  in  the 
future  phases  of  this  study.  The  details  of  this  TD-UTD  solution  will  be  submitted  for 
publication  soon  [21].  Some  numerical  results  are  shown  in  Figs.  21-22  for  a  source  at 
(x  =  a>  y  =  o,  z  =  0)  on  a  circular  cylinder  where  the  source  excited  by  a  transient  pulse 
which  is  the  same  as  in  Fig.  13.  The  circular  cylinder  is  chosen  because  an  exact  solution 
for  this  case  can  be  obtained  as  a  reference  solution  for  comparison  just  as  in  the  previous 
radiation  problem  case.  These  results  in  Figs.  21-22  are  obtained  via  an  ATT  based 
convolution  of  the  pulse  with  a  step  or  impulse  response,  and  the  TD-UTD  based  results 
are  compared  with  the  corresponding  exact  eigen  function  based  frequency  domain 
results  which  have  been  converted  into  the  time  domain.  The  results  indicate  the  surface 
fields  at  a  given  location  Q  (specified  by  the  angle  $).  The  agreement  between  the  TD- 
UTD  and  eigen  function  based  results  is  very  good. 
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4  Discussion 

A  set  of  five  FD-UTD  solutions  is  briefly  reviewed  in  this  report;  these  have  been 

* 

chosen  to  be  initially  incorporated  into  the  WE  algorithm  which  is  being  developed  by 
Monopole  Research.  The  subroutines  for  computing  the  FD-UTD  transition  functions  and 
hence  the  diffraction  functions  have  been  made  available  to  Monopole  Research.  The 
corresponding  five  TD-UTD  solutions  are  also  briefly  discussed.  These  five  solutions  in 
each,  the  FD-UTD,  and  also  the  corresponding  TD-UTD  representations,  involve  the 
diffraction  at  an  arbitrary  curved  wedge,  the  slope  diffraction  at  an  arbitrary  curved 
wedge,  the  surface  diffraction  by  a  smooth  convex  surface,  the  radiation  by  a  source  on  a 
smooth  convex  surface,  and  the  mutual  coupling  between  antennas  on  a  smooth  convex 
surface,  respectively.  It  is  noted  that  the  latter  two  constitute  new  solutions  in  the  TD- 
UTD  case  which  have  been  completed  under  this  contract.  These  FD-UTD  and 
corresponding  TD-UTD  solutions  when  implemented  within  the  WE  algorithm  can 
efficiently  analyze  the  radiation/scattering  /coupling  of  waves  in  complex  structures. 
Hence,  these  solutions  are  very  basic  and  useful  in  the  analysis  of  scattering,  radiation 
and  EMC/EMP  problems  of  practical  interest.  Additional  important  FD-UTD  results  are 
planned  to  be  included  within  the  WE  algorithm  in  the  future  phases  of  this  study  to  make 
the  WE  algorithm  more  versatile  in  applications  to  realistic  problems;  the  corresponding 
TD-UTD  solutions  will  also  be  developed  in  the  future  phases. 
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